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Although  it  is  accepted  that  the  significant  eccentricity  of  Mercury  (0.206)  favours  entrapment  into  the 
3:2  spin-orbit  resonance,  open  are  the  questions  of  how  and  when  the  capture  took  place.  A  recent  work 
by  Makarov  (Makarov,  V.V.  [2012],  Astrophys.  J„  752,  73)  has  proven  that  trapping  into  this  state  is  cer¬ 
tain  for  eccentricities  larger  than  0.2,  provided  we  use  a  realistic  tidal  model  based  on  the  Darwin-Kaula 
expansion  of  the  tidal  torque. 

While  in  Ibid,  a  Mercury-like  planet  had  its  eccentricity  fixed,  we  take  into  account  its  evolution.  To 
that  end,  a  family  of  possible  histories  of  the  eccentricity  is  generated,  based  on  synthetic  time  evolution 
consistent  with  the  expected  statistics  of  the  distribution  of  eccentricity.  We  employ  a  model  of  tidal 
friction,  which  takes  into  account  both  the  rheology  and  self-gravitation  of  the  planet. 

As  opposed  to  the  commonly  used  constant  time  lag  (CTL)  and  constant  phase  lag  (CPL)  models,  the 
physics-based  tidal  model  changes  dramatically  the  statistics  of  the  possible  final  spin  states.  First,  we 
discover  that  after  only  one  encounter  with  the  spin-orbit  3:2  resonance  this  resonance  becomes  the 
most  probable  end-state.  Second,  if  a  capture  into  this  (or  any  other)  resonance  takes  place,  the  capture 
becomes  final,  several  crossings  of  the  same  state  being  forbidden  by  our  model.  Third,  within  our  model 
the  trapping  of  Mercury  happens  much  faster  than  previously  believed:  for  most  histories,  10-20  Myr  are 
sufficient.  Fourth,  even  a  weak  laminar  friction  between  the  solid  mantle  and  a  molten  core  would  most 
likely  result  in  a  capture  in  the  2:1  or  even  higher  resonance,  which  is  confirmed  both  semi-analytically 
and  by  limited  numerical  simulations. 

So  the  principal  novelty  of  our  paper  is  that  the  3:2  end-state  is  more  ancient  than  the  same  end-state 
obtained  when  the  constant  time  lag  model  is  employed.  The  swift  capture  justifies  our  treatment  of  Mer¬ 
cury  as  a  homogeneous,  unstratified  body  whose  liquid  core  had  not  yet  formed  by  the  time  of  trapping. 

We  also  provide  a  critical  analysis  of  the  hypothesis  by  Wieczorek  et  al.  (Wieczorek,  M.A.,  Correia, 
A.C.M.,  Le  Feuvre,  M„  Laskar,  J.,  Rambaux,  N.  [2012],  Nat.  Geosci.,  5,  18-21)  that  the  early  Mercury  might 
had  been  retrograde,  whereafter  it  synchronised  its  spin  and  then  accelerated  it  to  the  3:2  resonance. 
Accurate  processing  of  the  available  data  on  cratering  does  not  support  that  hypothesis,  while  the 
employment  of  a  realistic  rheology  invalidates  a  key  element  of  the  hypothesis,  an  intermediate  pseudo- 
synchronous  state  needed  to  spin-up  to  the  3:2  resonance. 

©  2014  Elsevier  Inc.  All  rights  reserved. 


1 .  Motivation  and  plan 

Half  a  century  ago,  radar  observations  determined  the  mercuri- 
an  spin  period  to  be  «58  days  (Pettengill  and  Dyce,  1965),  which 
corresponds  to  a  3:2  spin-orbit  resonance.  A  later  study  (Margot 
et  al.,  2007)  revealed  that  the  orientation  of  Mercury’s  spin  axis 


*  Corresponding  author.  Fax:  +32  81724914. 

E-mail  addresses:  benoit.noyelles@unamur.be  (B.  Noyelles),  frouard@imcce.fr 
(J.  Frouard),  vvm@usno.navy.mil  (V.V.  Makarov),  michael. efroimsky@usno.navy.mil 
(M.  Efroimsky). 

http://dx.doi.Org/10.1016/j.icarus.2014.05.045 
0019-1035/©  2014  Elsevier  Inc.  All  rights  reserved. 


is  consistent  with  the  Cassini  State  1  (Colombo,  1965),  the 
obliquity  being  2.04  ±  0.08'  (Margot  et  al.,  2012). 

This  raised  the  question:  how  had  Mercury  been  trapped  into 
this  resonance?  A  consensus  exists  that  the  high  eccentricity  of 
Mercury  (currently  e  ss  0.206)  favours  the  trapping  by  widening 
the  resonance.  At  the  same  time,  there  is  a  cleavage  in  opinion 
on  whether  the  3:2  resonance  was  the  most  likely  end-state  of 
Mercury’s  spin-orbit  evolution. 

The  first  work  on  the  despinning  of  Mercury  was  published  by 
Goldreich  and  Peale  (1966).  They  obtained  a  7%  probability  for 
the  capture  into  the  3:2  resonance,  assuming  that  the  eccentricity 
of  Mercury  has  always  had  its  present  value,  and  that  the  tidal 
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torque  has  obeyed  the  MacDonald  (1964)  model.  Some  40  years 
later,  Correia  and  Laskar  (2004)  explored  a  chaotic  evolution  of 
Mercury’s  eccentricity,  showing  that  repetitive  episodes  of  eccen¬ 
tricity  increases  could  have  boosted  the  probability  of  entrapment 
in  the  3:2  resonance  to  55%. 

A  more  complete  study,  though,  should  include  a  large  fluid 
core.  Keeping  the  eccentricity  constant,  Peale  and  Boss  (1977a) 
explained  that  the  core-mantle  friction  could  have  significantly 
enhanced  the  3:2  capture  probability,  provided  that  Mercury  had 
not  been  trapped  into  the  2:1  resonance  prior  to  that.  Correia 
and  Laskar  (2009)  revisited  the  problem,  with  the  core  included, 
and  suggested  that  the  eccentricity  of  Mercury  might  have  been 
very  low  in  the  past,  lowering  drastically  the  probability  of  an  early 
capture  into  the  2:1  resonance,  and  making  a  later  trapping  into 
the  3:2  spin  state  almost  certain.  Both  works  by  Correia  and 
Laskar  (2004,  2009)  were  based  on  the  CTL  (constant  time  lag)  tidal 
model. 

Wieczorek  et  al.  (2012)  hypothesised  that  Mercury  had  initially 
been  retrograde,  evolving  later  to  synchronism.  Correia  and  Laskar 
(2012)  achieved  the  synchronism  from  an  initially  prograde  rota¬ 
tion.  Then  this  resonance  could  have  been  destabilised  by  a  huge 
impact  driving  Mercury  into  a  state  of  stable  pseudosynchronous 
rotation  maintained  by  a  significant  eccentricity,  before  falling  into 
the  3:2  resonance.  The  stability  of  pseudosynchronous  spin 
depends  on  the  applicability  of  the  CTL  model. 

Unfortunately,  neither  the  MacDonald  nor  the  CTL  models  are 
compatible  with  the  properties  of  telluric  bodies  (Efroimsky, 
2012a).  A  study  based  on  a  more  realistic  tidal  response  was  under¬ 
taken  recently  by  Makarov  (2012)  for  a  uniform  Mercury  analogue. 
The  study  demonstrated  that  entrapment  in  the  3:2  resonance  is 
inevitable  at  eccentricities  between  0.2  and  0.41,  without  invoking 
a  core-mantle  friction  or  eccentricity  variation. 

We  continue  this  line  of  research,  for  realistic  histories  of  the 
eccentricity.  In  Section  2,  we  overview  the  preceding  studies. 
Section  3  deals  with  the  triaxiality-generated  torque.  Section  4  pre¬ 
sents  an  expression  for  the  tidal  torque,  with  a  special  emphasis  on 
how  an  essentially  arbitrary  rheology  should  enter  that  expression. 
Section  5  addresses  the  despinning  timescale.  In  Section  6,  we 
compute  histories  of  Mercury’s  eccentricity.  Section  7  describes 
our  numerical  method.  We  test  it  by  reproducing  the  results  by 
Correia  and  Laskar  (2004). 

In  Section  8,  we  apply  our  method  to  a  homogeneous,  originally 
prograde  Mercury,  and  demonstrate  that  its  likeliest  end-state  is 
3:2.  For  a  cold  early  Mercury  (with  the  Maxwell  time 
tm  =  500  yr),  it  is  reached  in  less  than  20  Myr.  For  a  warmer 
Mercury  (with  lower  values  of  rM),  it  is  reached  in  less  than 
10  Myr.  This  justifies  the  homogeneity  assumption,  because  heat- 
ing-up  and  differentiation  takes  up  to  a  billion  of  years  (see  Appen¬ 
dix  A).  Despite  this,  to  make  our  study  complete,  in  Section  9  we 
explore  the  despinning  of  an  initially  prograde  Mercury  with  a 
liquid  core.  In  this  case,  even  a  moderate  friction  in  the  core-man¬ 
tle  boundary  boosts  the  probabilities  of  capture  into  the  3:2  and 
higher  resonances.  For  example,  with  the  present  value  of  Mer¬ 
cury's  eccentricity,  entrapment  into  the  2:1  resonance  becomes 
certain  (while  trapping  in  the  5:2  resonance  becomes  more 
probable  than  traversing  it).  Therefore,  had  Mercury  been  hot 
and  differentiated  soon  after  its  creation,  its  despinning  all  the 
way  to  the  3:2  state  would  be  highly  unlikely.  In  Section  10,  we 
provide  a  critical  analysis  of  the  hypothesis  by  Wieczorek  et  al. 
(2012)  that  the  early  Mercury  might  had  been  retrograde, 
whereafter  it  synchronised  and  then  accelerated  its  spin  to  the 
3:2  resonance.  Accurate  processing  of  the  available  data  on 


1  Besides  this,  the  MacDonald  model  is  inherently  inconsistent  (Williams  and 

Efroimsky,  2012;  Efroimsky  and  Makarov,  2013)  and  cannot  be  used  even  in 
theoretical  studies. 


cratering  does  not  support  that  hypothesis;  while  the  employment 
of  a  realistic  rheology  invalidates  a  key  element  of  that  hypothesis, 
an  intermediate  pseudosynchronous  state  needed  to  spin-up  to  the 
3:2  resonance.  So  that  hypothesis  has  to  be  abandoned.  Finally,  in 
Section  10  we  draw  conclusions. 

This  paper  is  an  extremely  squeezed  version  of  a  more  compre¬ 
hensive  work  (Noyelles  et  al.,  2013)  to  which  we  refer  here  and 
there  for  details  and  derivations. 

2.  Scenarios  of  Mercury’s  spin-orbit  evolution.  An  overview 

Consider  a  planet  of  a  mean  radius  R,  mass  Mptanet  and  the  prin¬ 
cipal  moments  of  inertia  A  <  B  <  C.  Assume  that  the  spin  of  the 
planet  is  directed  along  its  major-inertia  axis  z,  the  one  related  to 
the  maximal  moment  of  inertia  C.  The  sidereal  angle  0  of  the  planet 
can  be  reckoned  from  the  line  of  apsides  to  the  largest-elongation 
axis  x,  the  one  related  to  the  minimal  moment  of  inertia  A,  as  in 
Fig.  1.  The  star  exerts  two  torques  on  the  planet.  One,  T(7RI),  is 
due  to  the  planet’s  permanent  triaxiality.  Another,  T(nD£>,  is  caused 
by  the  tidal  deformation  of  the  planet.  Rotation  of  the  planet  about 
its  major-inertia  axis  z  is  then  governed  by  the  torques’  polar  com¬ 
ponents.  Expressing  the  maximal  moment  of  inertia  as 
C  =  £MpimetR1 2,  we  write  the  equation  of  motion  as 

-j-(TRI)  /j-(TIDE)  /j-(TRI)  ^-(TIDE) 

e=  C  iMplanetR2  (1) 

The  value  of  £  reflects  the  degree  of  inhomogeneity.  While  for  a 
homogeneous  sphere  { =  0.4,  for  Mercury  it  is  known  to  be 
0.346  ±0.014  (Margot  et  al.,  2012;  Noyelles  and  Lhotka,  2013). 
The  fact  that  £  <  0.4  indicates  that  the  inner  part  of  Mercury  is  den¬ 
ser  than  its  outer  part. 

This  equation  governs  the  spin  history.  The  outcome  also 
depends  on  the  initial  conditions  and  on  disruptive  events.  Three 
options  have  been  discussed  in  the  literature  hitherto. 

2.1.  Option  1:  Spin-down  of  a  homogeneous  Mercury 

Suppose  Mercury  was  captured  when  it  was  yet  homogeneous, 
and  its  liquid  core  was  formed  afterwards.  Simplistic,  this  approach 
is  convenient  mathematically,  wherefore  it  was  used  in  the  pioneer 
works.  The  first  such  investigation  was  undertaken  by  Goldreich 
and  Peale  (1966)  who  assumed  that  the  initially  spin  was  prograde 


Fig.  1.  The  principal  axes  x  and  y  of  the  planet  correspond  to  the  minimal  and 
middle  moments  of  inertia,  respectively.  The  horizontal  line  is  that  of  apsides,  so/  is 
the  true  anomaly.  In  neglect  of  the  apsidal  precession,  0  is  the  sidereal  angle  of  the 
planet.  The  angle  i (/  =f  -  0  renders  the  separation  between  the  planetocentric 
direction  towards  the  star  and  the  minimal-inertia  axis  x. 


28 


B.  Noyelles  et  al./ Icarus  241  (2014)  26-44 


and  pretty  fast.  Their  study  produced  a  7%  probability  of  capture 
into  3:2  resonance.  Therefore,  it  was  accepted  that  the  current 
rotational  state  of  Mercury  was  an  outcome  of  a  rather  unlikely 
event. 

Counselman  (1969)  proposed  that  the  eccentricity  of  Mercury 
could  have  varied  under  the  influence  of  other  planets,  and  that 
a  larger  value  could  enhance  the  probability  of  capture  into  the 
3:2  resonance.  However,  he  saw  no  reason  why  the  eccentricity 
would  have  been  bigger  than  0.25.  The  question  remained  open 
for  a  quarter  of  a  century,  until  Correia  and  Laskar  (2004)  ran 
1000  numerical  simulations  of  the  orbital  evolution  of  Mercury 
over  4  Gyr.  They  then  extracted  1000  histories  of  Mercury’s 
eccentricity  and,  among  other  things,  found  that  it  could  have 
reached  0.45.  The  authors  integrated  the  Eq.  (1 )  with  the  following 
expression  for  the  mean  tidal  torque  (sometimes  referred  to  as  the 
‘viscous  model’): 


D5 

r<nD£)  =  -3n2 *Mstark2At-^ \6A(e)  -  nAf(e)J , 


A(e)  = 


1  +  3e2  +§e4 
(1  -  e2)9/2 


1  -I-  i^P- 

W(e)  =  +  2 


15  p2  i  45  o4 


(i  -  e2) 


(2) 

(3) 


R,  n  and  e  being  Mercury’s  radius,  mean  motion  and  eccentricity; 
Mstar  being  the  mass  of  the  Sun;  the  Love  number  being  k2  =  0.4; 
and  the  time  lag  being  chosen  (pretty  arbitrarily)  as  nAf  =  1  /50. 
Integrations  began  with  a  spin  period  of  20  days,  to  see  how  many 
histories  lead  to  the  current  period  of  87.9  days,  i.e.,  to  the  3:2 
spin-orbit  resonance.2 

The  formulae  (2)  and  (3)  were  based  on  the  assumption  that  the 
time  lag  At  bears  no  dependence  upon  frequency.  Hence  the  name: 
the  constant  time  lag  (CTL)  model.  Detailed  derivation  of  these  for¬ 
mulae  can  be  found  in  the  Appendix  to  Williams  and  Efroimsky 
(2012).  In  various  forms,  this  result  appeared  much  earlier  in  Hut 
(1981)  and  Eggleton  et  al.  (1998).  In  an  implicit  form,  it  was 
pioneered  by  Goldreich  and  Peale  (1966,  Eq.  (24)). 

An  important  attribute  of  the  CTL  model  is  a  pseudosynchro- 
nous  equilibrium  spin  rate 

6  =  n  +  6ne2  +  ^ne4  +  0(e6).  (4) 

8 

For  the  present  eccentricity  of  «0.206,  this  would  be  9  rj  1.26n.  If 
the  eccentricity  reaches  0.285,  the  pseudosynchronous  equilibrium 
spin  rate  reaches  the  3:2  spin-orbit  resonance.  This  dramatically 
enhances  the  trapping,  as  the  planet  can  get  captured  and  can  stay 
in  this  resonance  even  after  the  eccentricity  decreases  to  a  lower 
value.  After  many  such  crossings,  the  overall  probability  of  capture 
into  the  3:2  state  must  exceed  1/2. 

The  physical  validity  of  this  scenario  is  questionable,  because 
the  existence  of  the  pseudosynchronous  spin  state  depends  on 
the  tidal  model  employed  -  and  because  the  CTL  model  is  very  dif¬ 
ferent  from  the  actual  behaviour  of  solids  and  partial  melts.  An 
accurate  description  of  tides  must  be  based  on  a  Fourier  decompo¬ 
sition  of  the  perturbing  potential,  with  each  term  corresponding  to 
an  appropriate  tidal  mode.  From  these  series,  a  similar  one  is 
derived  for  the  torque.  Then  the  tidal  response  at  each  mode 
should  be  entered.  Mathematically,  this  means  that,  into  each 
term,  an  appropriate  value  of  k2/Q  (or,  more  generally,  of  lq/Q,) 
must  be  inserted.  Insertion  of  realistic  frequency-dependencies  of 
ki/Qi  renders  pseudosynchronous  rotation  states  unstable,  as 
demonstrated  in  Makarov  and  Efroimsky  (2013). 


2  In  Correia  and  Laskar  (2004),  the  product  nAr  was  denoted  with  Q.  The  notation 

was  misleading,  as  the  so-defined  Q  did  not  coincide  with  the  quality  factor  (see 

Williams  and  Efroimsky,  2012,  Section  4). 


Basing  their  study  on  the  alleged  existence  of  the  pseudosyn¬ 
chronous  state  (4),  Correia  and  Laskar  (2004)  obtained  three  types 
of  capture: 

1.  Mercury  is  trapped  at  the  first  crossing  of  a  resonance  (31  tra¬ 
jectories  over  1000). 

2.  At  the  crossing  of  the  resonance,  Mercury  is  not  trapped,  but  the 
pseudosynchronous  spin  rate  (4)  turns  out  to  be  close  enough  to 
the  resonance,  to  induce  several  crossings  while  the  eccentricity 
oscillates;  and  the  planet  gets  locked  (168  trajectories). 

3.  Mercury  crosses  the  resonance  and  keeps  despinning  till  reach¬ 
ing  the  pseudosynchronous  rotation.  Then  secular  variations  of 
Mercury’s  eccentricity  respin  it,  because  increase  in  eccentricity 
entails  increase  in  the  pseudosynchronous  spin  rate.  The  planet 
gets  into  the  3:2  resonance  or  crosses  it  (to  be  despun  back 
later,  when  the  eccentricity  becomes  lower).  Multiple  crossings 
of  this  kind  can  result  in  a  trapping  (355  trajectories). 

22  trajectories  terminated  in  the  2:1  resonance,  554  in  the  3:2  res¬ 
onance,  and  36  in  the  synchronous  state.  As  the  current  3:2  reso¬ 
nance  appeared  to  be  the  most  probable  end  state,  the  problem 
of  Mercury’s  spin-orbit  resonance  seemed  to  have  been  solved. 
As  we  explained  above,  the  presented  solution  critically  depends 
on  the  choice  of  the  tidal  model. 

2.2.  Option  2:  Spin-down  of  a  differentiated  Mercury 

Based  on  the  magnetic  field  measurements  by  Mariner  10  (Ness 
et  al,  1974),  Mercury  was  suspected  to  harbour  a  partially  molten 
or  completely  liquid  core  (Ness  et  al.,  1975).  The  core  was  detected 
by  Earth-based  radar  observations  of  the  longitudinal  librations 
(Margot  et  al.,  2007),  which  are  significantly  larger  in  amplitude 
than  what  should  be  expected  from  a  uniformly  solid  planet.  As 
a  consequence,  extra  dissipation  should  be  taking  place  at  the 
core-mantle  boundary  (CMB).  If  the  core  had  formed  before  Mer¬ 
cury  was  trapped  into  the  3:2  resonance,  it  must  be  taken  into 
account  in  computation  of  the  despinning.  So  the  question 
becomes:  when  did  the  planet  acquire  its  core?  An  answer  to  this 
question  depends  upon  the  scenario  of  formation  of  Mercury.  Any 
such  theory  must  be  fit  to  explain  why  Mercury  is  poor  in  silicates, 
a  fact  ensuing  from  the  planet’s  high  density. 

One  scenario  is  that  of  early  volatilisation.  Within  this  approach, 
the  planet’s  composition  may  be  explained  if  the  removal  process 
in  Mercury’s  zone  of  the  solar  nebula  was  only  slightly  more  effec¬ 
tive  for  silicates  than  for  iron.  The  fractionation  required  to  pro¬ 
duce  an  iron-rich  planet  is  achieved  through  a  combination  of 
the  gravitational  and  drag  forces  acting  in  the  nebula  during  the 
formation  of  proto-Mercury  (Weidenschilling,  1978).  Another  pos¬ 
sibility  is  a  slow  volatilisation  scenario  which  suggests  that  after 
the  formation  of  Mercury  the  solar  wind  could  have  partially  volat¬ 
ilised  the  mantle  (Cameron,  1985;  Fegley  and  Cameron,  1987). 
Within  both  these  scenarios,  differentiation  of  the  planet  was  grad¬ 
ual.  As  explained  in  Appendix  A,  core  formation  required  hundreds 
of  millions  to  a  billion  of  years. 

A  different  scenario  is  rapid  differentiation  through  an  impact 
energetic  enough  to  have  molten  a  chondritic  protoplanet  (Benz 
et  al.,  1988, 2007).  The  event  would  have  resulted  in  a  hot  tumbling 
Mercury  with  a  CMB  forming  quickly.  The  hypothesis,  however,  is 
not  supported  by  the  MESSENGER  data.  As  was  shown  by 
Peplowski  et  al.  (2011),  the  abundance  of  potassium  is  inconsistent 
with  those  physical  models  of  Mercury  formation,  which  require 
extreme  heating  of  the  planet.  This  observation  prompted  (Wurm 
et  al.,  2013)  to  explore  the  action  of  photophoretic  forces  on  irradi¬ 
ated  solids.  The  study  revealed  that  silicates  are  preferentially 
pushed  into  an  optically  thick  disc.  The  subsequent  planetesimal 
production  at  the  outward-moving  edge  yields  metal-rich 
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planetesimals  close  to  the  star  and  metal-depleted  ones  farther  out 
in  the  nebula.  After  the  early  separation  of  metals  and  silicates,  the 
development  of  planet  goes  gradually,  with  a  slow  emergence  of  a 
core-mantle  boundary. 

Goldreich  and  Peale  (1967)  introduced  the  core-mantle  friction 
to  explain  the  despinning  of  Venus.  The  friction  was  modelled  as  a 
torque  linear  in  the  relative  angular  velocity.  Counselman  (1969) 
included  that  torque  into  the  theory  of  despinning  of  Mercury, 
and  concluded  that  the  core-mantle  friction  increased  the  local 
probability  of  capture  (the  probability  of  capture  into  each 
resonance  considered  separately).  The  new  torque  did  not  change 
considerably  the  global  probability,  in  that  Mercury  was  still  most 
likely  to  end  up  in  the  2:1  spin-orbit  resonance  and  had  very  few 
chances  to  reach  the  3:2  one.  This  result  was  confirmed  and  sup¬ 
plemented  by  Peale  and  Boss  (1977a)  who  concluded  that  the  tidal 
Q  should  be  smaller  than  100,  and  that  the  kinematic  viscosity  of 
the  molten  core  with  a  laminar  boundary  layer  should  be  compa¬ 
rable  to  that  of  water  (wO.Ol  cm2  s_1).  The  role  of  core-mantle  fric¬ 
tion  was  revisited  by  Correia  and  Laskar  (2009)  who  stated  that, 
due  to  the  chaotic  evolution  of  the  eccentricity  of  Mercury,  the 
planet  might  have  been  trapped  into  the  2:1  resonance  before  a 
decrease  of  the  eccentricity  disrupted  this  configuration.  In 
Correia  and  Laskar  (2009),  the  3:2  end-state  has  a  probability  of 
26%,  the  2:1  resonance  being  the  most  probable  outcome.  Anyway, 
the  authors  state  that  the  probability  of  a  final  capture  in  the  3:2 
state  can  be  increased  up  to  55%  or  73%  if  the  eccentricity  of  Mercury 
in  the  past  has  descended  below  the  critical  values  0.025  or  0.005, 
respectively.  Such  a  low  eccentricity  can  be  reached  through  con¬ 
sidering  short-period  effects  omitted  within  a  secular  model 
(Laskar,  2008).  All  of  these  studies  have  been  conducted  using 
the  CTL  tidal  model  which  permits  pseudosynchronism.  However, 
contrarily  to  the  rigid-Mercury  case,  this  time  the  pseudosynchro- 
nous  equilibrium  state  does  not  help  to  get  trapped  into  the  3:2 
resonance,  since  one  crossing  is  usually  enough. 

2.3.  Option  3:  Post-synchronous  evolution 

The  counts  of  craters  on  the  surface  of  Mercury,  from  the  Mar¬ 
iner  10  data  («45%  coverage)  and  the  first  two  flybys  of  MESSEN¬ 
GER,  led  Wieczorek  et  al.  (2012)  and  Correia  and  Laskar  (2012)  to 
suggest  that  Mercury’s  rotation  was  initially  retrograde,  and  the 
tides  (possibly  helped  by  the  core-mantle  friction)  accelerated  it 
to  the  1:1  resonance.3  While  in  the  synchronous  state,  Mercury 
acquired  much  of  its  craters,  whereafter  its  synchronism  was  desta¬ 
bilised  by  a  huge  impact.  The  fortuitously  directed  impact  drove 
Mercury  to  faster  spin  rates.  Assisted  by  a  high  eccentricity,  the  pla¬ 
net  reached  the  state  of  pseudosynchronous  rotation  which  is  stable 
in  the  framework  of  the  CTL  model.  The  chaotic  evolution  of  Mercury’s 
eccentricity  afterwards  got  it  finally  trapped  as  in  Correia  and  Laskar 
(2004). 

As  demonstated  in  Section  10,  an  impact  causing  the  disruption 
of  the  synchronous  state  should  have  left  a  crater  of  at  least  250- 
400  km  in  diameter,  depending  on  the  velocity  and  the  density  of 
the  impactor.  For  a  synchronous  Mercury  to  reach  directly  the 
3:2  resonance  without  help  from  tides  or  planetary  perturbations, 
the  crater’s  diameter  should  have  been  at  least  600  km. 

The  hypothesis  by  Wieczorek  et  al.  (2012)  is  based  on  three 
assumptions:  that  the  early  Mercury  was  impacted  heavily,  that 
the  crater  counts  are  evidence  of  a  past  synchronous  rotation,  and 
that  there  is  a  stable  state  of  pseudosynchronous  rotation.  Above 
that,  this  scenario  implies  that  Mercury  was  initially  retrograde. 


3  We  would  add  that,  as  an  alternative  possibility,  Mercury  might  had  started  in  a 
prograde  mode  and  got  decelerated  by  tides  all  the  way  down  to  synchronism, 
crossing  the  higher  resonances.  This  possibility,  however,  is  much  less  probable,  as  it 
requires  a  very  small  eccentricity  at  the  start  of  Mercury’s  history. 


Several  studies  confirm  that  the  early  Mercury  was  impacted 
intensively,  in  particular  during  the  Late  Heavy  Bombardment 
(Strom  et  al.,  2008,  2011).  The  associated  planetesimals  probably 
originated  from  the  Main  Asteroid  Belt  that  was  disturbed  by  the 
inward  migration  of  Jupiter  some  3.9-3.8  Ga  ago,  when  the  age 
of  Mercury  was  about  0.7  Ga  (Gomes  et  al.,  2005);  the  resulting 
craters  are  known  as  Population  1.  There  is  a  Population  2  of  youn¬ 
ger  craters  that  could  be  dominated  by  Yarkovsky-effect-driven 
near-Earth  asteroids  (Morbidelli  and  Vokrouhlicky,  2003;  Strom 
et  al.,  2005).  These  impacts  probably  triggered  an  episode  of  volca- 
nism  (Elkins-Tanton  and  Hager,  2005)  that  was  global  during  the 
LHB  and  occured  only  over  small  patches  or  within  large  impact 
basins  (Marchi  et  al.,  2013).  It  is  unclear  as  to  what  extent  this  voi- 
canism  played  a  role  in  internal  melt  generation.  From  the  global 
distribution  of  large  impact  basins  (Fassett  et  al.,  2012),  15  craters 
bigger  than  600  km  are  certain  or  probable,  while  21  are  suggested 
but  unverified.  If  we  set  the  limit  at  300  km  (collisions  strong 
enough  to  leave  the  synchronous  state  but  not  to  reach  the  3:2  res¬ 
onance  without  external  help),  then  we  have  46  certain  or  probable 
impact  basins,  and  61  suggested  but  unverified. 

The  counts  of  craters  in  Wieczorek  et  al.  (2012)  are  nearly  pre- 
MESSENGER,  in  that  they  combine  the  Mariner-10  data  with  those 
from  MESSENGER’S  first  two  flybys,  not  with  those  obtained  after 
the  orbital  insertion.  In  Fassett  et  al.  (2012),  the  inhomogeneous 
distribution  of  craters  is  emphasised,  and  the  authors  admit  that 
it  is  consistent  with  a  past  synchronous  spin.  They  also  mention 
another  possibility:  differential  resurfacing  of  Mercury  would 
explain  the  non-uniform  distribution  of  smooth  plains  (Denevi 
et  al.,  2009)  that  might  have  buried  degraded  basins.  Similarly, 
Strom  et  al.  (2011)  say:  "the  crater  size-frequency  distribution  on 
Mercury  shows  pronounced  regional  variations  that  we  argue  are 
primarily  the  consequence  of  differences  in  the  extent  and  age  of 
volcanic  plains  emplacement". 

The  initially  retrograde  spin  of  Mercury  is  an  interesting  option. 
Dones  and  Tremaine  (1993)  showed  that  collisions  affect  planets’ 
obliquities.  Kokubo  and  Ida  (2007)  got  a  distribution  for  their  val¬ 
ues.  Correia  and  Laskar  (2010)  pointed  that  an  oblique  Mercury  is 
less  likely  to  be  trapped  in  a  spin-orbit  resonance  than  a  non-obli- 
que  one.  So  most  studies,  including  ours,  still  start  with  the  obliq¬ 
uity  damped  to  0°  (prograde  spin)  or  180°  (retrograde  spin). 

Finally,  we  would  reiterate  that  the  stability  of  a  pseudosyn¬ 
chronous  equilibrium  depends  on  the  tidal  model  chosen.  For 
physics-based  rheologies,  pseudosynchronous  spin  is  transient 
and  cannot  be  used  to  help  the  angular  velocity  grow  from  the 
1:1  to  3:2  state.  This  gap  has  to  be  covered  in  one  jump,  which 
requires  a  significantly  larger  impactor. 


3.  The  triaxiality-caused  torque 

The  polar  torque  will  be  approximated  with  its  quadrupole  part 
(see,  e.g.,  Danby,  1962): 

T™  =l(B-A)  Sin  2^  - 1  (B  -  A)n2  sin  2(0  -/) .  (5) 

As  usual,  G  is  the  Newton  gravitational  constant,  Mstar  is  the  mass  of 
the  star,  r  is  the  distance  between  the  centres  of  mass,  while  yk  is  the 
angle  between  the  planetocentric  direction  towards  the  star  and  the 
minimal-inertia  axis  x,  as  in  Fig.  1.  This  angle  is  equal  to  the  differ¬ 
ence  between  the  planet’s  true  anomaly/,  as  seen  from  the  star,  and 
the  planet’s  sidereal  angle  9,  as  measured  from  the  line  of  apsides. 
These  and  other  notations  are  listed  in  Table  1,  the  numerical  values 
being  presented  in  Table  2.  The  mean  motion  is  defined  as  n  =  M, 
though  yGlMpia^f+M^YJa^  often  remains  a  good  approximation. 
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Table  1 

Symbol  key. 

Notation  Description 

f  The  moment  of  inertia  coefficient  of  Mercury 

R  The  radius  of  Mercury 

j-fTRi)  The  polar  component  of  the  triaxiality-caused  torque  acting  on 
Mercury 

'j-cnDE)  The  polar  component  of  the  tidal  torque  acting  on  Mercury 
Mpianet  The  mass  of  Mercury 

Mstar  The  mass  of  the  Sun 

a  The  semimajor  axis  of  Mercury 

r  The  instantaneous  distance  between  Mercury  and  the  Sun 

/  The  true  anomaly  of  Mercury 

e  The  orbital  eccentricity 

M  The  mean  anomaly  of  Mercury 

C  The  maximal  moment  of  inertia  of  Mercury 

B  The  middle  moment  of  inertia  of  Mercury 

A  The  minimal  moment  of  inertia  of  Mercury 

n  The  mean  motion 

G  Newton’s  gravitational  constant 

tm  The  viscoelastic  characteristic  time  (Maxwell  time) 

zA  The  inelastic  characteristic  time  (Andrade  time) 

/!  The  unrelaxed  rigidity  modulus 

J  The  unrelaxed  compliance 

a  The  Andrade  parameter 


Table  2 

The  numerical  values  of  the  parameters  entering  the  model.  The  value  of  the 
dimensionless  moment  of  inertia  is  taken  from  (Margot  et  al.,  2012),  and  the 
triaxiality  (B  -A)/C  is  derived  from  the  gravity  field  coefficient  C2 2  =  8.088  x  10-6 
(Smith  et  al.,  2012).  The  value  of  Mercury’s  radius  R,  was  borrowed  from  (Archinal 
et  al.,  2011).  The  values  of  both  the  semimajor  axis  a,  and  the  mean  motion  n,  are 
secular,  obtained  from  the  secular  planetary  theory  of  Bretagnon  (1982). 


Parameter 

Notation 

Numerical  value 

Semimajor  axis 

a 

5.791  x  107  km 

Mean  motion 

n 

26.0879  rad/yr 

Radius  of  Mercury 

R 

2440  km 

Dimensionless  moment  of  inertia 

(  =  C/(Mpl0„etR2) 

0.346 

Triaxiality  of  Mercury 

(B-A)/C  =  4C22/£ 

9.350  x  10~5 

Ratio  of  masses  of  the  Sun  and 
Mercury 

M star /M planet 

6.0276  x  106 

Unrelaxed  rigidity 

p 

8  x  1010  Pa 

Maxwell  time 

tM 

500  yr 

Andrade  time 

Ta 

500  yr 

Andrade  parameter 

a 

0.2 

Newton’s  gravitational  constant 

G 

66468  m3  kg-'  yr2 

To  spare  ourselves  of  the  necessity  to  compute  the  true  anomaly 
/,  we  shall  express  the  torque  through  the  mean  anomaly  M.  To  do 
so,  we  shall  rewrite  the  expression  (5)  as 

3  a3 

T™  =  -  -  (B  -  A)n2  (sin  26  sin  2/  -  cos  26  cos  2/)  (6) 

and  shall  insert  therein  the  following  formulae  (see,  e.g.,  Duriez, 
2002): 

/A\ 3  +°° 

(-)  cos(2/)  =  X^3,2(e)  cos(kM),  (7) 

k=— 00 

Q  sin(2/)  =  ^Xk3'2(e)sin(fcX),  (8) 

k=— oo 

Xk,m(e)  being  the  classical  Hansen  coefficients.  These  are  computed 
through  the  formula 

1  °°  P 

K'm(e)  =  (1  +zTn~  E(-z)PEc^m+,c! lm+Jk-m+r-2h(l<e),  (9) 

p=0  h= 0 


where  z=  (1  -  v/l  -  e2)/e,  while  Cab  are  the  binomial  coefficients, 
and  Jk(x)  are  the  Bessel  function  of  the  first  kind.  The  Hansen 
coefficients  are  related  to  the  eccentricity  functions  via 
G/p,  =  X,l%y-2p.  Altogether,  formulae  (6)-(8)  entail: 

=  _|(B_A)n 2  g  C20,(e)sin  (2[fl-  (l  +  §)wi]).  (10) 

q=—oo 

In  practice,  the  infinite  series  may  be  truncated  to  a  sum  over 

q  =  -4, _ 6.  To  our  study,  most  relevant  resonances  are  the  1:1, 

3:2,  2:1,  5:2,  3:1,  7:2  and  4:1  ones,  i.e.,  those  corresponding  to 

q  =  0 . 6  for  the  prograde  rotation,  and  the  resonances  -2:1, 

-3:2,  -1:1,  -1:2,  1:2  and  1:1  for  the  retrograde  rotation,  i.e. 

q  =  -6, _ 0.  In  the  latter  case,  the  terms  corresponding  to 

q  =  - 6,-5  have  been  added  and  proved  to  be  unsignificant 
(Section  10).  For  e  >  0.3,  computation  of  the  eccentricity  function 
Gipq(e)  requires  a  lot  of  computer  time,  as  the  convergence  of  the 
series  (9)  slows  down.  So  we  computed  these  functions  once  for 
the  whole  range  of  eccentricities,  and  interpolated  them  during 
the  computation,  using  cubic  splines.  To  that  end,  the  GNU 
Scientific  Library  (Galassi  et  al.,  2009)  was  employed. 

4.  The  tidal  torque 

Tidal  torques  acting  on  a  telluric  body  have  long  been  subject  to 
misconceptions.  One  of  them  has  it  that  certain  rheologies  should 
be  regarded  “unphysical".  The  grounds  for  this  belief  were  that 
T<tide>  a  k2/Q>  wherefore  a  quality  factor  Q  scaling  as  a  positive 
power  of  the  tidal  frequency  would  render  an  infinite  torque  in 
the  zero-frequency  limit.  The  oversight  stemmed  from  confusing 
the  seismic  quality  factor  Q  with  the  tidal  quality  factors  defined 
as  1  /Q,  =  sin  e(,  where  e(  are  phase  lags.  The  seismic  Q  factor 
reflects  the  rheology,  while  the  tidal  Q  reflects  the  interplay  of  rhe¬ 
ology  and  self-gravitation,  a  circumstance  that  excludes  unphysical 
divergencies  in  the  quality  functions  Jq/Qj  =  /q(tO/mpg)  sinei(mimp?). 

A  torque  obtained  in  the  CTL  model,  was  employed  by  some 
authors  in  the  constant  angular  lag  context.  A  study  of  this  error 
can  be  found  in  Efroimsky  and  Makarov  (2013). 

Another  common  misconception  is  the  belief  in  “pseudosyn- 
chronous”  spin  of  terrestrial  planets  and  moons.  While  the  stability 
of  such  states  for  bodies  with  oceans  remains  an  open  question,  for 
telluric  objects  such  states  are  unstable  and  therefore  transient,  as 
explained  in  Makarov  and  Efroimsky  (2013).  The  explanation  can 
be  extended  also  to  molten  planets,  except  in  the  case  of  extremely 
low  viscosity. 

4.1.  General  facts 

A  consistent  way  of  treating  linear  tides  comprises  two  steps: 

(1 )  Both  the  tide-raising  potential  of  the  secondary  body  (in  our 
case,  of  the  star)  and  the  incremental  tidal  potential  of  the 
distorted  primary  (in  our  case,  of  the  planet)  must  be 
expanded  into  Fourier  series.  The  terms  in  both  series  are 
numbered  by  four  integers  Impq. 

(2)  An  Impq  term  in  the  latter  series  is  related  to  the  Impq  term 
in  the  former  series,  via  the  dynamical  Love  number  kt  and 
phase  lag  £/  as  functions  of  the  Fourier  mode  coimpq. 

An  Impq  term  of  the  induced  tidal  potential  is  proportional 
to  the  Impq  term  of  the  tide-raising  potential,  multiplied 
by  some  geometric  factor  and  by  the  function 
k/cosq  =  fc|(ftj|mp,)cose/(ct)lmp,).  Similarly,  an  Impq  term  of 
the  secular  part  of  the  tidal  torque  will  be  proportional  to 
the  Impq  term  of  the  tide-raising  potential,  multiplied  by  a 
geometric  factor  and  by  the  so-called  quality  function 
ki  sin  £i  =  ki(a>impq)  sin  ei(a>impq).  These  results,  now  classical, 
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were  obtained  by  Kaula  (1964).  For  a  significant  precursor, 
see  Darwin  (1879).  As  was  pointed  out  later  by  Efroimsky 
(2012a),  the  torque  also  contains  an  oscillating  part. 

4.2.  An  expression  for  the  tidal  torque 

A  Fourier  expansion  of  the  tidal  torque  was  written  down  by 
Goldreich  and  Peale  (1966)  who  took  into  account  only  its  secular 
part.  A  schematic  proof  of  that  formula  was  presented  later  by 
Dobrovolskis  (2007).  A  comprehensive  derivation  of  the  polar  com¬ 
ponent  of  the  torque  can  be  found  in  Efroimsky  (2012a)  where  it  is 
demonstrated  that  the  torque  contains  both  a  secular  and  a  rapidly 
oscillating  part. 

An  infinite  series,  that  expression  can  be  truncated  for  a 
planet  which  is  not  too  close  to  the  star  (R/a  <c  1),  has  a  small 
obliquity  (i  ~  0),  and  not  very  large  an  eccentricity  (e  <  0.4). 
All  three  conditions  are  obeyed  by  Mercury  now,  and  we  have 
no  reason  to  suspect  that  they  were  ever  violated  in  the  past.4 
The  appropriate  truncation  is  carried  out  step-by-step  in  the 
Appendix  B  to  Noyelles  et  al.  (2013),  and  it  renders  the  following 
expression: 

rfM)  =  3  CMLfxyf  C2oq(e)G2oj(e)k2( m220,) 

k  /  qj=- 1 

x  sine2(ciJ220ij)  cos  [(q  - j)M ]  +  0(e8e)  +  0(i2e) 

+  0((,R/a)7e),  (11a) 

where  R  and  0  are  the  radius  and  rotation  angle  of  the  planet,  0  is  its 
spin  rate,  a  and  e  are  the  semimajor  axis  and  eccentricity,  Gipq(e)  are 
the  eccentricity  polynomials. 

An  Impqhj  term  of  the  torque  contains  a  quality  function 
ki(coimpq)  sin  ei(coimpq).  The  dynamical  Love  number  l<i(coimpq)  is  a 
positive  definite,  even  function  of  the  mode.5  The  phase  lag 
£i (coimPq)  is  an  odd  function  whose  sign  always  coincides  with 
that  of  the  tidal  mode  a>impq  (e.g.,  Efroimsky  and  Makarov,  2013). 
Thus  each  factor  ki(a>impq)  sin  ei(coimpq)  can  be  written  as  ki(a>impq) 
sin  | e/ ( co(mpq )  | Sgn  (a> impq).  The  terms  with  q  ¥=  j  are  oscillating.  While 
they  influence  the  fate  of  each  trajectory,  the  overall  statistics  is 
defined  overwhelmingly  by  the  secular  part  (Makarov  et  al.,  2012). 
The  terms  with  q  =  j  constitute  the  secular  part: 

(7'TIDE)>,  2=^GMlrR5a-6J2G220q(e)k2(co220q) 

q-  i 

x  sin  |e2(m220ii)|  Sgn  (m220q)  +  0(e8e) 

+  0(i2e).  (lib) 

In  neglect  of  the  apsidal  and  nodal  precessions,  the  Fourier  tidal 
modes  (atmpq  look  as 

coimpq  =  (l-2p  +  q)n-md,  (12) 

where  l>  2  is  an  integer  termed  the  degree,  m  =  0 . /  is  an 

integer  known  as  the  order,  while  p  and  q  are  two  other  integers, 
p  varying  from  0  through  I,  and  q  taking  values  from  -oo  to  oo. 
Each  quality  function  can  be  expressed  as  a  function  of  the 
planetary  spin  0: 


4  In  our  modelling  of  Mercury’s  orbital  evolution,  the  eccentricity  reached  the 
values  between  e  =  0.4  and  e  =  0.45  in  only  four  histories  out  of  a  thousand,  and  over 
very  limited  time  spans.  This  happened  for  three  simulations  covering  0.5  Gyr  and  for 
one  simulation  covering  0.1  Gyr.  The  obliquity  of  Mercury  never  exceeded  3.3  arcsec 

(Peale,  2005). 

5  So  it  would  be  equally  legitimate  to  write  ki(a>imm)  as  where 

Ximpq  —  l(0imp?l  are  the  positive  definite  physical  forcing  frequencies  corresponding 
to  the  Fourier  modes  oj/mp,. 


ki(o)inipq)  sin  \ti((j)jmpq)  Sgn(co/mp,) 

=  ki((l  -2p  +  q)n  -  mb)  sin  |ef((Z  -  2p  +  q)n 
-  mi)) |  Sgn((l  -  2p  +  q)n  -  mi)),  (13a) 

where  ((/  -  2p  +  q)n  -  mb)  denotes  a  functional  dependence  on  the 
argument  (I  -  2p  +  q)n  -  mb,  not  multiplication  by  a  factor  of 
(/  -  2p  +  q)n  -  mb.  In  our  case,  (13a)  becomes: 

k2(to22Qq)  sin  |e2(co220g)|  Sgn  (co22oq) 

=  k2(2(n  -  9)  +  qn)  sin  |e2(2(n  -  0)  +  qn)  |  Sgn(2(n  -  9) 

+  qn).  (13b) 

In  each  term  of  the  sums  (11),  the  factor  k2sine2  = 
k2  sin  |e2|Sgn((u220p),  expressed  as  a  function  of  6,  has  the  shape  of 
a  kink,  as  shown  in  Fig.  2.  The  shape  ensues  from  interplay  of  rhe¬ 
ology  (which  dominates  at  higher  frequencies)  and  self-gravitation 
(which  takes  over  in  the  zero-frequency  limit),  see  the  Sections  4.3 
and  4.4.  The  emergence  of  this  shape  is  natural,  for  each  term 
should  transcend  zero  and  change  its  sign  continuously  when  the 
rotation  rate  goes  through  the  appropriate  spin-orbit  resonance.6 

In  the  sum  (11),  the  kink-shaped  factors  k2  sine2  are  multiplied 
by  G20p(e).  So  the  sum,  treated  as  a  function  of  0,  is  a  superposition 
of  nine  kinks  of  different  magnitudes.  The  kinks  are  centred  at  dif¬ 
ferent  resonant  values  of  0.  The  resulting  curve  depicting  the  sum 
(11)  crosses  zero  in  points  very  close  to  the  resonances 
6  =  n(l  +q/2),  but  not  exactly  in  the  resonances  —  like  the  little 
kink  near  9/n  =  3/2  in  Fig.  3.7 

Sitting  on  the  right  slope  of  a  more  powerful  kink  (the  one  with 
q  =  0),  the  little  kink  (corresponding  to  q  =  1)  experiences  a  verti¬ 
cal  displacement  called  bias.  Although  the  bias  of  the  little  kink  is 
caused  by  all  the  other  kinks,  the  influence  of  its  left  neighbour 
is  clearly  dominant.  The  vertical  bias  of  the  little  kink  in  Fig.  3  is 
a  feature  immanent  to  any  spin-orbit  resonance,  if  the  eccentricity 
is  non-zero.  The  zero-eccentricity  case  is  exceptional  because,  in 
the  quadrupole  approximation,8  it  comprises  only  the  semidiurnal 
component  of  the  tidal  torque  (the  one  corresponding  to  the 
principal  tidal  frequency  oj220o). 

A  constant  bias  of  the  secular  torque  was  also  considered  by 
Goldreich  and  Peale  (1968),  who  pointed  out  its  crucial  role  in 
the  capture  probability  calculation.  However,  in  our  model  the  bias 
has  a  physical  meaning  very  different  from  that  in  the  CTL  model. 
In  our  case,  the  bias  is  generated  by  the  inputs  from  the  non¬ 
resonant  modes  called  into  being  by  non-uniform  orbital  motion. 
Without  the  non-vanishing  bias  terms,  capture  into  a  spin-orbit 
resonance  would  be  inevitable,  and  the  tidal  dissipation  would 
completely  cease  once  the  rotator  is  locked  in  the  resonance.  Near 
a  resonance  q',  the  spin  rate  0  is  close  to  n(l  +  q'/2),  and  the  right- 
hand  side  of  the  sum  (lib)  can  be  conveniently  decomposed  into 
two  parts.  The  first  one  is  the  q  =  q1  term,  a  kink-shaped  function 
transcending  zero  at  0  =  n(l  +  q1/ 2).  The  second  part,  bias, 
comprises  all  the  q  ^  q'  terms  of  the  sum.  This  way,  the  bias  is 
the  input  of  all  the  q  ^  q'  terms  into  the  values  assumed  by  the 
overall  torque  in  the  vicinity  of  the  q  =  q'  resonance.  In  the  said 


6  In  a  hypothetical  situation  of  a  planet  despinning  through  a  resonance  at  a 
constant  rate,  the  appropriate  Fourier  mode  is  linear  in  time.  Then  the  tidal  torque, 
expressed  as  a  function  of  time,  will  acquire  a  similar  kink  shape.  This  situation  is 
described  in  Ferraz-Mello  (2013,  Fig.  7b).  Any  physically  acceptable  rheological  model 
must  lead  to  this  or  similar  kind  of  tidal  torque  behaviour  near  a  resonance. 

7  As  ensues  from  (11a)  and  (13a),  the  Impq  term  of  the  secular  ( j  =  q)  part  of  the 
torque  is  decelerating  for  ro/mm  <  0  and  accelerating  for  wimpq  >  0.  Equivalently,  it  is 
negative  on  the  right  of  the  Impq  resonance  and  positive  on  its  left  -  see  the  two  kinks 
depicted  in  Fig.  3. 

8  From  the  tables  of  the  F(i)  and  C(e)  functions,  it  can  be  observed  that  only  the 
terms  with  q  =  p  =  l  —  m  —  0  come  out  non-vanishing  for  e  =  i  =  0.  For  1  =  2,  the  sole 
surviving  term  is  {Impq}  =  {2200},  while  for  /  >  2  we  would  get  the  extra  modes  with 
{Impq}  =  {1100},  for  each  1  included. 
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Fig.  2.  Atypical  shape  of  the  quality  function  k:(ca)  sinei(<o),  where  co  is  a  shortened 
notation  for  the  tidal  Fourier  mode  cotepq. 

vicinity,  the  bias  is  a  slowly  changing  function,  which  can,  to  a 
good  approximation,  be  treated  as  constant.  For  not  too  large 
eccentricities,  the  bias  is  usually  negative  -  this  is  why  the  little 
kink  in  Fig.  3  is  shifted  downwards. 

While  the  q  =  q'  term  assumes  a  zero  value  at  6/n  =  (1  +  q'/2), 
the  bias  displaces  the  location  of  the  zero.  The  3:2  kink  in  Fig.  3 
goes  through  nil  in  a  point  located  a  bit  to  the  left  of  0/n  =  3/2. 
Being  extremely  small,  the  shift  of  the  equilibrium  away  from 
the  resonance  has  no  practical  consequences,  since  the  net  nonzero 
tidal  torque  is  compensated  by  an  opposing  secular  triaxiality- 
caused  torque,  as  explained  in  Makarov  and  Efroimsky  (2013). 

4.3.  The  tidal  torque  and  the  role  of  rheology 

The  torque  (11),  taken  as  a  function  of  the  spin  8,  is  a  superpo¬ 
sition  of  kinks.  Depicting  a  typical  kink  and  the  way  kinks  superim¬ 
pose,  Figs.  2  and  3  serve  illustrative  purposes  and  do  not 
correspond  to  realistic  parameters.  Each  kink  is  stretched  horizon¬ 
tally,  its  extrema  being  spread  apart,  to  emphasise  that  the  transi¬ 
tion  through  nil  is  continuous. 

Analogous  figures  for  realistic  parameters  of  Mercury  would 
have  sharper  features,  like  in  Fig.  4.  Despite  the  sharp-looking 
peaks,  the  dependency  is  continuous,  and  the  middle  part  has  a 
finite  slope,  even  though  in  the  picture  it  looks  almost  vertical. 


Fig.  3.  The  angular  acceleration  due  to  the  secular  part  of  the  tidal  torque. 
Expanded  into  series,  the  torque  is  expressed  as  a  function  of  the  dimensionless 
spin  rate,  B/n.  The  left  kink  comes  from  the  leading,  semidiurnal  term 
{Impq}  =  {2200}.  It  is  centred  in  the  point  B/n  =  1  corresponding  to  the  1:1 
spin-orbit  resonance.  The  right  kink  comes  from  the  term  with  {Impq}  =  {2201}, 
which  is  vanishing  in  the  3:2  resonance.  The  right  kink  is  considerably  displaced 
vertically,  owing  to  the  bias  provided  by  the  slopes  of  other  kinks,  mainly  by  the 
leading,  semidiurnal  kink. 


rate  of  rotation  9/ n 

Fig.  4.  Angular  acceleration  of  Mercury  caused  by  the  secular  tidal  torque  (11),  in 
the  vicinity  of  the  2:1  spin-orbit  resonance. 


The  form  of  the  plot  is  dictated  by  the  shape  and  magnitudes  of 
the  quality  functions  ki(oJimpq)  sin  ei(a>impq)  expressed  as  functions 
of  the  spin  rate  via  formulae  (13).  The  quality  functions  show  up 
in  each  term  of  the  tidal  torque,  and  the  shapes  of  these  functions 
are  crucial  for  determining  the  probabilities  of  entrapment.  The 
shapes  are  defined  by  the  size  and  mass  of  the  body  and  by  its  rhe¬ 
ological  properties.  The  rheological  properties  are  accounted  for  by 
a  constitutive  equation,  i.e.,  by  a  law  interconnecting  the  strain  and 
stress.  Within  linear  rheologies,  a  constitutive  equation  can  be 
rewritten  in  the  frequency  domain  where  each  Fourier  harmonic 
uyv(yj  of  the  strain  gets  expressed,  algebraically,  as  the  appropriate 
harmonic  of  the  stress  ffyv(x)  multiplied  by  the  complex 
compliance  J(%)  appropriate  to  this  harmonic: 

2  Uyv(x)  =J(X)Vyv(X),  04) 

X  =  Ximpq  being  the  physical  forcing  frequency,  and  yv  being  the  ten- 
sorial  indices.  The  relation  (14)  is  an  exact  analogue  to  a  complex 
harmonic  of  current  written  as  the  appropriate  harmonic  of  voltage, 
multiplied  by  the  inverse  complex  impedance.  The  roles  of  the  cur¬ 
rent,  voltage  and  inverse  impedance  are  played  by  2uyv(x),  <7,V(Z) 
and  J(x),  correspondingly. 

The  complex  compliance  ](x)  carries  all  information  about  the 
response  of  the  material,  under  the  assumption  that  the  material 
is  isotropic.  Modelling  the  mantle  rheology,  i.e.,  deriving  ](x),  we 
rely  on  a  combined  model  which  is  Andrade  at  higher  frequencies 
and  Maxwell  in  the  low-frequency  limit.  As  was  demonstrated  in 
Efroimsky  (201 2a, b),  the  mathematical  formalism  of  the  Andrade 
model  can  be  modified  in  a  manner  permitting  a  switch  to  the 
Maxwell  model  at  low  frequencies.  The  necessity  for  this  switch 
originates  from  the  fact  that  at  different  frequencies  different 
physical  mechanisms  dominate  the  internal  friction. 

4.4.  ki/Qi  as  functionals  of  the  complex  compliance  of  the  mantle 

The  expression  (lib)  contains  the  so-called  quality  function 
ki((Oimpq)  sin  e,(oj,mp,)  =k,(co,mpq)  sin  \e,(coimpq)\Sgn(wlmpq)  defined  by 
the  rheology  and  self-gravitation  of  the  planet.  As  shown  in 
Efroimsky  (201 2a, b),  the  absolute  value  of  this  function  is 
calculated  as 

3  Ail 

ki(a>imPq )  sin  \ti(oJimpq)\  =  ~  (yj  +  a  )2  +  j2  ’  (^5) 

where  the  dimensionless  real  and  imaginary  parts  of  the  compli¬ 
ance  ](x)  of  the  mantle  are 

U  =  nizejix)]  =  1  +  (rr„r“cos  (f)r(a+l),  (16) 

1  =  plm\J(x)}  =  -(XTm)-1  -  (xta)-*  sin  (^)r( oh-  1);  (17) 


B.  Noyelles  et  al./ Icarus  241  (2014)  26-44 


33 


r  is  the  Gamma  function,  while  %  is  a  shortened  notation  for  the 
physical  forcing  frequency 

Ximpq  =  \(Oimpq\  =  |(l  -2 p  +  q)n  -  m0\,  (18) 

a  is  the  Andrade  parameter  (about  0.2  for  silicates  with  partial 
melt),  while  zM  and  zA  are  the  Maxwell  time  and  the  Andrade  time 
of  the  mantle.  The  quantities  A(  are  given  by 

A  3  fi(2l2 +41  +  3)  3(2f2  +41  +  3) 

4lnGp2R 2  4/7 iGp2R2J  ' 

where  p,  J  =  1  / p,  p  and  R  are  the  unrelaxed  rigidity,  compliance, 
mean  density  and  radius  of  the  planet,  while  G  is  the  gravitational 
constant.  The  presence  of  the  factors  (19)  in  the  denominator  of 
(15)  is  owed  to  the  competition  between  rheology  and  self¬ 
gravitation  (the  latter  playing  effectively  the  role  of  extra  rigidity). 

The  Maxwell  time  is  the  ratio  of  viscosity  to  rigidity: 
tm  =  y\l p  =  t] j.  The  Andrade  time  zA  was  introduced  in  Efroimsky 
(2012a, b)  as  a  substitute  to  the  thitherto  used  Andrade  parameter 
p  which  had  fractional  dimensions  and  thus  lacked  a  clear  physical 
interpretation.  Referring  the  reader  to  Efroimsky  (201 2a, b)  for 
details,  we  would  only  recall  here  that  beneath  some  threshold  fre¬ 
quency  the  mantle’s  response  becomes  close  to  that  of  a  Maxwell 
body.  Mathematically,  this  means  that  below  the  threshold  the 
parameter  zA  increases  rapidly  as  the  frequency  goes  down.  As  a 
result  of  this,  only  the  first  terms  in  (16)  and  (17)  stay  relevant. 

At  frequencies  near  and  above  the  threshold,  both  the 
viscoelastic  and  inelastic  mechanisms  are  present.  This  implies 
that  the  parameters  zM  and  zA  should  have  comparable  values,  as 
can  be  seen  from  the  above  expressions  (16)  and  (17)  for  the  com¬ 
pliance.  In  that  frequency  band,  the  mantle  behaves  as  an  Andrade 
body.  With  increase  of  the  frequency,  the  inelastic  reaction  takes 
over  the  viscoelastic  reaction,  as  defect-unpinning  becomes  a 
dominating  mechanism  of  friction. 

We  set  zA  =  tm  at  the  frequencies  above  the  threshold  (chosen 
to  be  1  yr-1,  like  for  solid  Earth).  Below  the  threshold,  we  set  zA  to 
increase  exponentially  with  the  decrease  of  the  frequency,  so  at 
low  frequencies  the  rheology  approaches  the  Maxwell  one.  Numer¬ 
ical  runs  show  that  the  ensuing  capture  probabilities  are  not  very 
sensitive  to  how  fast  the  switch  from  the  Andrade  to  Maxwell 
model  is  (Makarov,  2012). 

To  avoid  numerical  complications  near  the  zero  frequency  (i.e., 
on  crossing  resonances),  we  multiplied  both  the  numerator  and 
denominator  of  (15)  by  when  performing  integration. 

Numerical  implementation  of  the  tidal  model  (11)  takes  a  much 
longer  computation  time  than  the  model  (2).  So  we  run  a  set  of 
1000  simulations,  with  the  tidal  torque  containing  only  the  secular 
terms  (with  q  =  j).  With  various  Maxwell  times,  we  do  this  for  a 
rigid  Mercury  initially  prograde  (Section  8)  and  retrograde 
(Section  10).  For  a  differentiated  Mercury  with  core-mantle  fric¬ 
tion,  we  present  a  semi-analytical  model  combined  with  limited 
numerical  tests  -  for  the  planet  initially  prograde  (Section  9)  and 
retrograde  (Section  10). 

5.  A  characteristic  time  of  despinning 

Before  integrating  the  law  of  motion  (1),  it  is  of  both  theoretical 
and  practical  interest  to  estimate  the  characteristic  time  scale  of 
secular  deceleration.  This  time  scale  is  greatly  influenced  by  the 
tidal  dissipation  rate  which,  in  its  turn,  is  dependent  on  the  instan¬ 
taneous  rate  of  rotation,  on  the  ratio  of  planet’s  radius  R  to  the  orbi¬ 
tal  separation  a,  on  the  eccentricity  e,  and  on  other  physical 
parameters.  Theory  tells  us  that  some  bodies  spin  down  quickly, 
on  the  time  scale  of  105— 106  yr  (e.g.,  the  Moon),  while  others  need 
time  intervals  comparable  to  their  lifetime,  to  slow  down 
appreciably  (Castillo-Rogez  et  al.,  2011).  So  the  required 


Fig.  5.  Decimal  logarithm  of  Mercury's  characteristic  spin-down  time  as  a  function 
of  the  spin  rate,  for  e  =  0.1  (upper  curve),  e  =  0.2056  (middle  curve),  e  =  0.3  (lower 
curve).  The  curves  are  computed  for  discrete  sets  of  points,  to  avoid  the  resonances. 
For  all  values  of  e,  the  spin  is  accelerating  when  i)/n  <  1  (not  shown  in  this  figure). 
For  a  sufficiently  large  e,  the  spin  may  be  accelerating  also  over  some  spans  of 
values  of  8/n  larger  than  unity.  E.g.,  for  e  =  0.3  it  accelerates  when  8/n  e  [1, 1.5), 
wherefore  the  segment  of  Log[-8/0]  between  8/n  =  1  and  8/n  =  1.5  is  not  plotted 
here.  At  this  eccentricity,  the  secular  acceleration  8  is  positive  for  8/n  <  1.5. 

computing  time  may  be  prohibitively  long  for  a  definite  result  to 
be  obtained. 

Traditionally  (e.g.,  Rasio  et  al.,  1996),  a  characteristic  time  of 
despinning  is  defined  as 

tdespin  =  TTT7  •  (20) 

w 

Representing  a  time  span  over  which  the  spin  slows  considerably, 
this  quantity  should  not  be  taken  for  the  actual  time  a  planet  needs 
to  slow  down  from  a  given  0  to,  say,  the  synchronous  state.  The 
actual  time  of  despinning  can  be  obtained  only  by  numerical  inte¬ 
gration  of  (1). 

The  decimal  logarithms  of  characteristic  spin-down  times  are 
shown  in  Fig.  5  for  three  values  of  eccentricity:  e  =  0.1  (upper 
curve),  e  =  0.2056  (middle  curve),  and  e  =  0.3  (lower  curve).  Gen¬ 
erally,  these  times  are  of  the  order  of  107  yr  or  longer,  while  the 
integration  step  size  has  limitations  imposed  by  the  kink  shape 
of  the  tidal  torque  components.  ’ 

6.  Numerical  modelling  of  the  history  of  Mercury’s  eccentricity 

6.1.  Background 

Here  we  present  a  method  to  obtain  the  eccentricity  history  of 
Mercury  over  Gyr  timescales.  It  is  based  on  statistical  formulae 
developed  in  Laskar  (2008). 

The  chaotic  behaviour  of  the  Solar  System  was  demonstrated, 
with  averaging,  by  Laskar  (1989)  and,  with  no  averaging  involved, 
by  Sussman  and  Wisdom  (1992).  No  precise  solution  can  be 
obtained  by  integration  back  in  time  over  ~50  Myr  or  longer 

(Laskar  et  al.,  2011). 

Within  a  statistical  approach,  Laskar  (2008)  studied  the  possible 
orbital  history  of  Mercury,  using  averaged  planetary  equations.  His 


9  When  the  planet  traverses  a  low-order  resonance,  a  considerable  slow-down 
occurs  due  to  the  left  shoulder  of  the  appropriate  tidal  kink.  In  Figs.  3  and  4,  the  right 
shoulder  of  each  kink  is,  typically,  negative  (decelerating),  while  the  left  shoulder  of 
each  kink  is,  typically,  positive  (accelerating).  We  prefer  to  say  ‘typically’,  because 
there  may  exist  situations  where  a  kink  corresponding  to  a  higher  resonance  is 
located,  fully  or  partially,  below  the  horizontal  axis  -  see  the  small  kink  in  Fig.  3.  Even 
though,  despinning  leftwards  through  this  kink  will  result  in  a  decrease  of  the  overall 
despinning  action  of  the  tidal  torque. 
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Fig.  6.  Left:  the  Rice  distribution  function  for  T=  1  Gyr  from  Eq.  (21)  (line)  and  its  approximation  obtained  from  the  distribution  of  500,000  synthetic  orbits  evolved  for  1  Gyr 
with  the  Wiener  process  of  Eq.  (23)  (crosses).  Right:  evolution  over  time  of  five  synthetic  orbits. 


dynamical  model  contains  the  eight  planets,  with  their  mutual  sec¬ 
ular  gravitational  perturbations  obtained  through  a  second-order 
perturbation  method,  along  with  perturbations  from  the  General 
Relativity  and  the  Moon  (Laskar,  1986).  Note  that  the  averaged 
equations  imply  constant  secular  semi-major  axes  of  the  planets. 
After  numerical  integration  of  a  large  number  of  histories  of  the 
Solar  System  over  5  Gyr,  starting  with  very  close  initial  conditions, 
some  statistics  can  be  obtained  about  the  values  of  the  eccentricity 
and  inclination  of  each  planet  over  this  time  span.  In  particular,  the 
eccentricity  of  Mercury  e  statistically  follows  a  Rice  density 
function  (see  Laskar,  2008): 


/r,s(e)  =  ^2  exP 


(-4^) 


(21) 


where  the  parameter  a  obeys  the  approximation 

a2(T)  =  b0  +  bJ  (22) 

corresponding  to  a  diffusion  process,  with  b0  =  2.07  x  10~3, 
bi  =  1.043  x  1 0  3,  s  =  0.1875,  and  with  T  being  the  time  in  Gyr.  /„ 
is  the  modified  Bessel  function  of  the  first  kind. 

Mind  that  s  and  a  are  not  the  mean  and  standard  deviation  of 
the  Rice  distribution. 


6.2.  Method 


Our  goal  here  is  to  produce  a  synthetic  time  series  of  Mercury’s 
eccentricity,  based  on  the  density  function  (21).  We  use  a  simple 
Wiener  process  (Brownian  motion)  of  the  type 

e(T+  St)  =  e(T)  +  a(St)Se,  (23 ) 

where  e(0)  is  the  initial  eccentricity,  and  a  is  the  standard  deviation. 
The  values  of  e  are  separated  in  time  by  a  constant  step  dt,  while  5e 
are  independent  and  stationary  increments.  For  a  Wiener  process, 
the  increments  &e  obey  a  Gaussian  distribution  Af(0, 1). 

To  obtain  the  standard  deviation  needed  in  (23),  we  fit  a 
Gaussian  distribution 


/(e;e,<72) 


1 

\j2no2 


exp 


(e  -  e)2 
2a2 


(24) 


to  the  Rice  density  function  (21),  e  being  the  mean  value,  and  the 
standard  deviation  obeying  a2(t)  =  0.0009t.  We  then  applied  a 
few  corrections  to  the  time  series. 

First,  the  Rice  density  has  a  characteristic  drift  of  its  mean  value 
and  most  probable  value  over  time.  For  the  Gaussian  distribution 
to  still  maintain  a  good  fit  in  the  complete  time  interval,  this  is 


compensated  by  adding  a  small  drift  D  =  0.0023T  to  the  value  of 
the  eccentricity  at  a  time  T.  By  fitting  the  drift,  we  obtain  an  initial 
mean  e  =  0.1876,  to  be  used  also  for  the  initial  value  e(0). 

Secondly,  we  change  the  timescale:  T*  =  T  +  2.33  Gyr.  This  step 
is  necessitated  by  the  fact  that  the  first  250  Myr  of  integration 
computed  by  Laskar  (2008)  have  a  very  different  slope  of  <r2(t), 
and  have  not  been  taken  into  account  in  his  fit  of  (21)  -  see 
Laskar  (2008,  Fig.  15b).  To  be  able  to  fit  the  Rice  distribution  with 
its  current  parameters  for  all  values  of  T,  a  Gaussian  distribution 
with  a  standard  deviation  linear  in  time  has  to  start  at  an  initial 
time  of  -2.33  Gyr.  Also,  due  to  the  lack  of  information  about  the 
first  0.5  Gyr  of  evolution,  we  keep  only  the  part  corresponding  to 
the  interval  T*  e  [0.5,4]  Gyr. 

One  of  the  discrepancies  between  the  Rice  and  Gaussian  distri¬ 
butions  lies  in  the  distribution  of  very  low  eccentricity  values.  The 
Rice  density  is  defined  in  the  range  [0,  +oo],  and  its  distribution  for 
low  e  tends  to  be  more  linear  with  increasing  T.  The  Gaussian 
distribution  is  defined  on  [— oo,  +oo],  allowing  the  possibility  of 
negative  eccentricities,  and  is  symmetric  with  respect  to  its  mean 
value.  We  can  see  from  Laskar  (2008)  that  the  eccentricity  distribu¬ 
tion  of  Mercury  does  not  show  a  strong  linear  behaviour  for  low  e 
(which  is  not  the  case  for  the  other  planets  of  the  Solar  System). 
The  orbits  which  have  e  <  0  at  any  time  during  the  3.5  Gyr  time 
span  have  been  discarded  from  the  final  results,  with  the  side  effect 
of  slightly  depleting  the  distribution  of  low  e,  compared  to  a  non- 
modified  Gaussian  distribution.  0  We  have  checked  that  it  does  not 
introduce  any  bias,  since  this  reduced  distribution  matches  the  Rice 
distribution  for  low  e  even  more  closely.  Fig.  6  shows  the  Rice 
distribution  function  for  T  =  1  Gyr  computed  with  the  Eq.  (21),  com¬ 
pared  with  the  distribution  of  500,000  synthetic  orbits  obtained 
through  the  process  (23).  Finally,  as  the  secular  planetary  equations 
used  by  Laskar  (2008)  are  time-reversible,  we  reverse  the  time  direc¬ 
tion  of  the  process  (see  the  examples  of  synthetic  orbits  in  Fig.  6). 


7.  Numerical  simulation  of  Mercury’s  spin  rate 

7.1.  Methodology 

Our  next  goal  is  to  study  the  influence  of  the  realistic  tidal 
model  on  the  probability  of  Mercury’s  capture  into  spin-orbit 
resonances.  Makarov  (2012)  has  recently  accomplished  this  work 
for  a  constant  eccentricity.  Here  we  use  our  synthetic  trajectories, 


10  Synthetic  orbits  obtained  via  the  Brownian-motion  method  can  have  negative 
eccentricities.  While  the  Rice  distribution  prohibits  this,  the  Gaussian  distribution 
does  not. 
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to  include  the  statistical  variations  of  Mercury’s  eccentricity  in  the 
study.  We  will  build  statistically  relevant  samples  of  simulations  of 
the  despinning  of  Mercury,  with  two  kinds  of  tidal  model:  the  one 
from  Correia  and  Laskar  (2004)  and  our  model  (lib). 

We  integrate  numerically  the  Eq.  (1)  for  the  spin  rate  of  Mer¬ 
cury.  The  torques  entering  this  equation  depend  on  the  eccentricity 
e  which  is  calculated,  at  each  step,  from  the  synthetic  distribution. 
The  value  of  e  at  a  given  time  is  estimated,  by  cubic  spline  interpo¬ 
lation,  from  the  synthetic  evolution  trajectory  sampled  at  a  regular 
grid  of  time  points  (Dahlquist  and  Bjorck,  2008).  The  interpolation 
error  is  of  the  4th  order  (Sonneveld,  1969),  so  a  decrease  of  the 
sampling  step  of  the  trajectory  by  a  factor  of  a  reduces  the  interpo¬ 
lation  error  by  a  factor  of  rs  a4.  For  generating  eccentricities,  we 
used  the  interpolation  routines  from  the  GNU  Scientific  Library 
(Galassi  et  al.,  2009).  The  numerical  integrator  was  the  Adams- 
Bashforth-Moulton  lOth-order  predictor-corrector  scheme 
(Hairer  et  al.,  1987).  Through  the  simulations,  the  value  of  the 
mean  motion  was  kept  constant:  n  =  26.0878  rad/yr. 


7.2.  Test:  reproducing  the  results  from  Correia  and  Laskar  (2004) 

To  test  our  method,  we  reproduced  the  results  from  Correia  and 
Laskar  (2004),  by  running  10  sets  of  1000  histories  each,  all  starting 
with  a  spin  period  of  20  days.  As  stated  in  Correia  and  Laskar 
(2004),  the  choice  of  the  initial  spin  rate  is  not  critical,  and  high- 
order  resonances  are  unimportant.  To  save  the  computer  time, 
we  integrated  over  only  3  Myr,  which  is  long  enough  to  slow  down 
Mercury  and  to  reach  the  Type  1  captures.  Each  set  of  1000  simula¬ 
tions  used  the  same  set  of  1000  synthetic  eccentricity  trajectories, 
the  differences  between  our  10  sets  being  only  in  the  initial  value 
of  the  sidereal  angle  8.  We  use  these  10  sets  to  check  how  the  num¬ 
ber  of  trapped  trajectories  can  change  from  one  set  to  another, 
with  the  same  eccentricities.  We  used  the  same  values  of  the  phys¬ 
ical  parameters  as  in  Correia  and  Laskar  (2004).  Details  of  these 
simulations  are  described  comprehensively  in  the  extended 
version  of  this  paper,  Noyelles  et  al.  (2013). 

The  obtained  statistics  varied  significantly  among  our  10  sets,  as 
was  expected  due  to  the  Poisson  distribution  of  the  number  of 
positive  outcomes.  We  got  between  29  and  56  Type  1  captures  into 
the  3:2  state,  while  Correia  and  Laskar  (2004)  had  31.  So  our  results 
are  consistent.  Our  method  is  thus  valid  and  can  be  combined  with 
a  different  tidal  model. 


8.  Revisiting  the  Option  1:  a  rigid  Mercury,  initially  prograde 

Here  we  present  the  statistics  of  capture  into  spin-orbit  reso¬ 
nances,  based  on  the  realistic  tidal  model,  i.e.,  on  the  formula 
(lib)  with  the  tidal  response  (15)-(17)  plugged  in.  We  integrated 
the  law  of  motion  (1 ),  with  the  torques  (10)  and  (lib)  inserted,  and 
with  the  eccentricity  evolving  along  the  synthetic  trajectories 
described  in  Section  6  and  tested  in  Section  7.  In  this  section, 
Mercury  is  initially  set  prograde,  with  a  spin  period  of  20  days. 

First,  we  set  the  values  of  the  Maxwell  and  Andrade  times  of  the 
mantle,  tm  and  ta,  emerging  in  the  formulae  (15)— (17).  Our  choice 
of  values  for  the  Maxwell  time  of  Mercury  is  in  agreement  with  the 
range  proposed  by  Padovan  et  al.  (2014).11 

As  was  explained  in  the  Section  4.4,  their  values  should  be  com¬ 
parable.  We  set  them  equal,  with  the  important  caveat  that  ta 
increases  rapidly  below  some  threshold  frequency,  to  make  the 
low-frequency  reaction  purely  Maxwell. 


11  For  the  Maxwell  time,  Padovan  et  al.  (2014)  used  a  fixed  reference  value,  to  which 
the  actual  tm  was  related  by  some  rescaling  factor  dependent  on  the  temperature  and 
on  the  grain  size.  With  that  factor  taken  into  account,  the  actual  Maxwell  time  of 
Mercury  would  assume  values  from  several  months  to  about  1 50  years. 


For  the  Earth,  tm  «  500  yr.  However,  only  relatively  cold 
mantles  can  have  such  a  large  rM.  For  warmer  mantles,  the  values 
of  dozens  of  years  or  even  years  should  be  expected.  Indeed,  the  vis¬ 
cosity  i ]  depends  on  the  temperature  T  through  the  Arrhenius  law 
17  oc  exp (A*/RT),  where  the  gas  constant  is  R  =  8.3  J/(mol  K)  and 
the  activation  energy  for  silicate  rocks  is  A*  rs  6  x  105  J  mol-1.  At 
the  same  time,  the  rigidity  p  depends  on  T  slower,  until  T 
approaches  about  three  quarters  of  the  melting  temperature,  the 
latter  being  increased  by  the  confining  pressure. 


8.1.  First  simulations  with  zM  =  tA  =  500  yr 

We  present  the  results  of  1000  numerical  runs.  The  capture  sta¬ 
tistics  are  as  follows: 


•  resonance  4:1 

•  resonance  7:2 

•  resonance  3:1 

•  resonance  5:2 

•  resonance  2:1 

•  resonance  3:2 

•  resonance  1:1 


5  captures  (0.5%), 

21  captures  (2.1%), 
40  captures  (4%), 

103  captures  (10.3%), 
279  captures  (27.9%), 
444  captures  (44.4%), 

104  captures  (10.4%). 


The  4  remaining  histories  are  between  the  1:1  and  3:2  reso¬ 
nances  and  are  bound  to  be  trapped  into  the  former,  as  the  secular 
tidal  torque  is  negative,  i.e.,  decelerating.  The  histogram  of  these 
histories  is  presented  in  Fig.  7. 

We  conclude  that  the  3:2  resonance  is  the  most  probable  for 
the  Type  I  captures.  Other  types  of  captures  are  much  less  likely, 
as  the  tidal  despinning  drives  the  spin  rate  to  the  synchronous 
rotation  0  =  n,  making  more  than  one  crossing  of  the  resonance 
impossible.  As  the  evolution  of  the  eccentricity  will  not 
change  these  results  of  resonant  trapping,  we  can  safely  avoid 
the  necessity  to  numerically  integrate  the  spin  evolution  further 
in  time. 

Trapping  into  the  3:2  resonance  usually  requires  no  more  than 
20  Myr,  often  much  less. 


8.2.  Different  values  of  the  Maxwell  time  tm 

Mercury  might  have  been  warmer  in  its  early  life,  especially  if 
impacted.  This  would  imply  shorter  values  for  the  Maxwell  time. 
Approximating  the  tidal  torque  with  its  secular  part,  we  carried 
out  the  simulation  for  zM  =  5  yr  and  zM  =  15yr.  In  both  cases, 
we  set  the  Andrade  time  equal  to  the  Maxwell  time.  The  results 
are  presented  in  Table  3  and  Fig.  8. 

Shorter  Maxwell  times  favour  higher-order  resonances, 
especially  the  2:1  one.  As  the  actual  final  state  of  Mercury  is 
the  3:2  resonance,  a  longer  initial  Maxwell  time  (hundreds  of 
years)  is  a  likelier  option.  At  the  same  time,  we  should  keep  in 
mind  that  the  final  states  are  probabilistic.  With  small  but  not 
negligible  probability,  the  planet  can,  in  principle,  reach  the 
3:2  spin-orbit  resonance  even  with  an  initial  Maxwell  time 
smaller  than  20  years. 

Another  consequence  of  a  shorter  initial  Maxwell  time  is  an 
even  faster  tidal  despinning.  While  our  former  simulations 
with  tm  =  500  yr  reached  the  3:2  resonance  in  less  than 
20  Myr,  with  the  smaller  values  we  arrived  at  a  final  state  in 
less  than  10  Myr. 


12  A  change  in  the  temperature  varies  the  viscosity  and  the  Maxwell  time: 

AfAR.  Consider  a  planet  with  a  silicate  mantle  of  an  average  T  =  2300K 
and  tm  =  500  yr.  A  decrease  in  the  viscosity  and  the  Maxwell  time  by  9/10  (down  to 
tm  =  50  yr)  would  correspond  to  an  increase  of  the  temperature  by  AT  =  66  K. 
Heating  up  the  planet  by  another  hundred  degrees,  we  can  lower  its  tm  down  to 
several  years. 
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CTL 


Fig.  7.  Distribution  of  the  final  spin  rates  of  10,000  trajectories  integrated  using  the  CTL  model  (left)  and  the  model  (15)-(17)  (right),  with  zM  =  tA  =  500  yr.  For  the  CTL,  we 
considered  only  one  crossing  of  the  3:2  resonance.  (After  that,  most  trajectories  go  into  pseudosynchronous  rotation.)  For  our  model  (15)-(17),  the  3:2  resonance  is  the 
likeliest  end-state  and  is  reached  within  20  Myr. 


Table  3 

Entrapment  statistics  for  shorter  values  of  the  Maxwell  time. 


rM  =  5  yr 

tM  =  15  yr 

Resonance  4:1 

13  captures  (1.3%) 

7  captures  (0.7%) 

Resonance  7:2 

43  captures  (4.3%) 

22  captures  (2.2%) 

Resonance  3:1 

118  captures  (11.8%) 

104  captures  (10.4%) 

Resonance  5:2 

240  captures  (24%) 

177  captures  (17.7%) 

Resonance  2:1 

361  captures  (36.1%) 

368  captures  (36.8%) 

Resonance  3:2 

205  captures  (20.5%) 

284  captures  (28.4%) 

Resonance  1 : 1 

20  captures  (2%) 

38  captures  (3.8%) 

9.  Revisiting  the  Option  2:  Mercury  with  a  liquid  core,  initially 
prograde 

As  indicated  by  the  large  amplitude  of  longitudinal  forced  libra- 
tions  (Margot  et  al.,  2007),  Mercury  has  a  massive  molten  core 
decoupled  from  the  mantle.  In  Appendix  A,  we  explain  that  gradual 
formation  of  a  core  could  have  taken  up  to  billion  of  years.  The 
short  spin-down  times  obtained  in  this  paper  (10-20  Myr)  tell  us 
that  the  core  formed  after  the  capture  into  the  3:2  resonance.  There 
exist,  however,  scenarios  wherein  the  core  emerged  at  the  earliest 


stages  (e.g.,  after  a  collision).  Hence  the  need  to  consider  the  case 
where  the  capture  was  influenced  by  the  core.  That  the  friction 
in  the  boundary  layer  can  enhance  capture  was  stated  by  Correia 
and  Laskar  (2009)  and,  earlier,  by  Peale  and  Boss  (1977a). 

In  Peale  and  Boss  (1977a),  the  authors  argue  that,  in  the  laminar 
boundary  approximation,  capture  of  Mercury  into  the  2:1  state 
becomes  certain  even  for  weak  coupling.  On  examining  the  treat¬ 
ment  in  Peale  and  Boss  (1977a),  we  found  that  the  said  study  was 
based  on  a  “constant-Q”  model,  i.e.,  on  an  assumption  that  the  sec¬ 
ular  tidal  torque  is  frequency-independent.  Below  we  shall  demon¬ 
strate  that  the  calculation  of  capture  probabilities  in  the  presence  of 
a  core,  developed  by  Goldreich  and  Peale  (1966, 1967, 1968),  can  be 
generalised  to  any  realistic  frequency-dependent  tidal  model.  We 
shall  then  apply  this  generalisation  to  our  model  (lib). 

9.1.  The  theory  of  Goldreich  and  Peale  (1966,  1967,  1968) 

The  first  steps  of  our  analysis  follow  the  seminal  study  by 
Goldreich  and  Peale  (1966, 1967, 1968)  on  the  spin  state  of  Venus. 
The  theory  describes  the  behaviour  of  the  quantity 


Fig.  8.  Capture  statistics  with  the  secular  tidal  torque  (lib),  for  a  warmer  Mercury,  i.e.  with  a  shorter  Maxwell  time  than  before. The  main  difference  is  a  higher  probability  for 
the  2:1  resonance. 
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y  =  8  - 


(25) 


which  is  the  angle  of  libration  about  the  resonance  number  q' .  It  is  a 
resonance  where  0/n  =  1  +  q' / 2,  with  M  being  the  mean  motion, 
and  with  0  and  0  being  the  sidereal  angle  and  spin  rate  of  Mercury. 
Our  definition  (25)  of  the  libration  angle  is  consistent  with  that 
introduced  by  Goldreich  and  Peale  (1967,  1968).  It  differs  by  a  fac¬ 
tor  of  1  / 2  from  the  quantity  y  =  20  -  (2  +  q')M  employed  in 
Makarov  (2012)  and  Makarov  et  al.  (2012). 

For  a  laminar  boundary  between  the  core  and  the  mantle,  the 
equations  of  motion  are 

(T^>)  k 

ym  =  -tog  sin  2ym  +  ^-i — L-  (ym-yc),  (26a) 

'-m  l-m 

7c  =  (^  (7m  -7c),  (26b) 

where  (r<T,DE|)  is  the  secular  part  of  the  polar  tidal  torque  acting  on 
the  mantle,  k  is  a  friction  constant,  while  the  indices  m  and  c  desig¬ 
nate  the  mantle  and  core,  respectively  (so  Cm  and  Cc  are  the  largest 
moments  of  inertia  of  the  mantle  and  the  core).  The  quantity 
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1  1/2 


(27) 


is  the  frequency  which  small-amplitude  (sin2ym  rj  2ym)  librations 
would  have,  were  there  no  friction.  As  shown  in  Appendix  B,  the 
separatrix  between  rotation  and  libration  is  given  by  the  equation  ! 


ym  =  V2<a>o  cosym. 


(28) 


The  Eqs.  (26)  describe  a  non-harmonic,  driven  pendulum  with  fric¬ 
tion.  When  coupling  is  strong  (i.e.,  when  k  is  large),  the  core  is  par¬ 
ticipating  in  the  librational  oscillations.  Opposite  is  a  situation  with 
weak  friction,  which  is  the  case  for  Mercury;  in  that  situation  the 
core  does  not  have  time  to  react  to  the  high-frequency  changes  in 
the  rotation  of  the  mantle.  As  the  planet  approaches  the  q'  reso¬ 
nance  from  the  domain  of  positive  ym,  the  mantle  decelerates  much 
faster  than  the  core,  increasing  the  negative  difference  in  their  spin 
rates.  In  the  works  by  Goldreich  and  Peale  (1966, 1967, 1968),  a  cor¬ 
nerstone  assumption  was  that  the  core  is  in  a  dynamic  equilibrium 
with  the  time-average  tidal  torque  (which  was  implicitly  set  con¬ 
stant).  In  our  treatment,  we  shall  depart  from  this  assumption 
and  shall  generalise  the  treatment  to  an  arbitrary  form  of  the  fre¬ 
quency-dependence  of  the  torque. 


9.2.  Generalisation  of  the  Goldreich-Peale  theory 

In  the  expression  for  the  degree-2  secular  torque  (lib),  a  220q 
term,  when  expressed  as  a  function  of  the  Fourier  mode  a>2 2oq,  is 
antisymmetric.  As  Goldreich  and  Peale  (1966,  1967,  1968)  were 
using  simplistic  models  with  a  constant  time  lag  or  a  constant  Q, 
an  impq  term  in  their  theory  was  either  constant  or  linear  in  the 
tidal  mode.  In  realistic  models,  each  term  is  a  kink,  as  in  Fig.  2. 
When  the  formula  (12)  is  utilised  to  express  each  term  as  a 
function  of  d/n,  the  overall  sum  becomes  a  superposition  of  kinks 
of  different  sizes,  see  Fig.  3. 


13  Our  Eq.  (28)  coincides  with  a  formula  from  Goldreich  and  Peale  (1968),  but  is 
slightly  different  from  a  similar  equation  in  Makarov  (2012)  and  Makarov  et  al. 
(2012).  A  difference  in  the  numerical  factor  is  due  to  the  fact  that  the  y  defined  in 
Goldreich  and  Peale  (1968)  and  in  the  current  paper  differs  by  a  factor  of  1/2  from  the 
y  introduced  in  Makarov  (2012)  and  Makarov  et  al.  (2012)  -  see  the  paragraph  after 
the  formula  (25)  above.  In  our  equations,  we  also  keep  the  mass  factor  omitted  in 
Goldreich  and  Peale  (1968).  We  denote  the  frequency  (27)  with  </2w0  to  conform  to 
the  notation  from  Peale  and  Boss  (1977a,b). 


9.2.1.  The  tidal  torque  near  a  resonance 

In  the  vicinity  of  a  resonance  q',  the  secular  tidal  torque  can  be 
written  down  as 

(Tz'DE))i_2  =  W  +  V.  (29) 

Here  W  is  the  q  =  q'  term  of  the  sum  (lib): 

W  =  lcM2starR5a-6G220q,(e)k2(co220q,) 

x  sin  |e2(a>220q')|Sgn(cq220q')-  (30) 

A  kink  centred  around  the  point  6/n  =  1  +  q'/2,  this  term  is  a  func¬ 
tion  of  the  Fourier  mode  co220q G  which  is  the  negative  twice  of  ym: 

®220ij'  =  — 2ym.  (31) 

So  we  may  treat  W  as  a  function  of  ym: 

W  =  VV(ym)  =  -KG20qlk2(ym)  sin  |e2(7m)|Sgn(ym),  (32) 

K  being  a  positive  constant.  The  term  V  in  (29)  is  comprised  by  the 
rest  of  the  sum  (lib): 

V  =  K^^G20qk2(co22oq)  sin  |e2(co22Oij)|Sgn(co22o,])  +  0(e8e) 

q*q' 

+  0(i2e).  (33a) 

In  the  vicinity  of  the  q  =  q'  resonance,  the  tidal  modes  are  approx¬ 
imated  with  ft>220ij  =  0/22011'  +  (q  -  q')n  ss  (q  -  q')n,  whence  V  is 
approximated  as 

V  «  Kj^C20qk2((q  -  q')ri)  sin  |e2((q  -  q')n)|Sgn(<j  -  q').  (33b) 

q*q' 

In  the  said  vicinity,  V  is  a  slowly-changing,  nearly  constant  bias  gen¬ 
erated  by  the  non-resonant  terms.  This  can  be  understood  from 
Fig.  3  where  a  smaller  kink  on  the  right  is  residing  on  the  smooth 
slope  produced  by  other  kinks,  mainly  by  the  larger  kink  on  the  left. 

The  idea  of  the  decomposition  (29)  belongs  to  Goldreich  and 
Peale  (1966,  1967,  1968)  who  employed  it  for  calculation  of 
entrapment  probabilities,  within  the  two  simplest  tidal  models, 
that  of  a  frequency-independent  Q  and  that  of  a  frequency-inde¬ 
pendent  time  lag.  In  the  works  by  Makarov  (2012)  and  Makarov 
et  al.  (2012),  their  theory  was  generalised  to  an  arbitrary  fre¬ 
quency-dependence  of  the  tidal  torque.  That  generalisation,  how¬ 
ever,  was  carried  out  for  a  homogeneous  rotator.  So  our  goal  now 
is  to  embrace  the  case  with  a  core. 

9.2.2.  Dynamics  of  the  mantle  near  a  resonance 

The  mode-dependent  part  of  the  torque  changes  very  quickly  as 
the  mode  coimpq  goes  through  a  resonance.  For  weak  coupling,  the 
core  does  not  react  appreciably  to  sudden  accelerations  of  the 
mantle.  So  the  term  k(ym  -  yc)/Cm  in  the  Eq.  (26a)  is  small.  Numer¬ 
ical  runs  of  (26)  demonstrate  that  the  resonance  is  traversed 
quickly  (within  years). 

Near  a  resonance,  the  steady  part  of  ym  vanishes:  (ym)  =  0 
(while  the  oscillating  part  survives  due  to  free  librations).  At  the 
same  time,  the  steady  part  of  ym  does  not  vanish  -  on  the  contrary, 
it  becomes  large.  Time-averaging  of  the  Eqs.  (26a)  and  (26b)  yields: 

(7m)  =^~  +  ^~  7c  (34a) 

'-m  *-m 

7c  =  -^7c-  (34b) 

We  also  assume  that  near  resonances  the  mantle  is  in  a  dynamical 
equilibrium  with  the  constant  part  V  of  the  torque.14  So  the  friction 


14  The  core  is  not  in  equilibrium,  because  of  the  velocity  offset  and,  thus,  a  nearly 
constant  friction  force. 


38 


B.  Noyelles  et  al./ Icarus  241  (2014)  26-44 


force  stays  approximately  constant  during  the  short  passage,  and  the 
quantities  (yra)  and  yc  are  approximately  equal.  Hence, 

k  .  V  k  .  ■  ,  ,  ■  V  Cc 

-ryz=7^  +  7^yc  or,  equivalently  :  yc  =  -  —  — — — — . 

(35) 


Insertion  of  the  Eqs.  (28),  (29)  and  (35)  into  the  equation  of  motion 
(26a)  furnishes  an  equation  equivalent  to  the  Eq.  (27)  in  Peale  and 
Boss  (1977b): 


ym  =  -®o  sin  2ym 


Wf s  V 

Cm  Cm  +  Cc 


V2I<OJ0 

Cm 


cos  y 


m' 


(36) 


9.2.3.  The  probability  of  entrapment 

To  estimate  capture  probabilities,  it  is  sufficient  to  consider  two 
librations  bracketing  the  point  of  resonance  y  =  0.  As  explained  in 

Goldreich  and  Peale  (1968),  the  probability  is 

SE  being  the  total  change  of  kinetic  energy  at  the  end  of  the 
libration  below  the  resonance,  and  A£  being  the  offset  from  zero 
at  the  beginning  of  the  last  libration  above  the  resonance.  To 
obtain  A£,  the  torques  acting  on  the  mantle  should  be  integrated 
over  a  cycle  of  libration.  To  obtain  SE,  the  torques  should  be  inte¬ 
grated  over  two  librations  symmetric  around  the  resonance 
to 220q'  =  0.  As  a  result,  the  odd  part  of  the  tidal  torque  at  q  =  q' 
doubles  in  the  integration  for  SE,  whereas  the  bias  vanishes.  As 
can  be  understood  from  the  Eq.  (36),  the  core-mantle  friction 
has  two  manifestations  in  the  expression  for  the  secular  torque. 
The  even  (driving)  part  V  gets  attenuated  by  a  factor  of 
Cm/(Cm  +  Cc).  The  odd  (restoring)  part  W  is  boosted  by  an  addi¬ 
tional  friction  term  proportional  to  kco0.  Both  these  effects  lead 
to  higher  probabilities  of  entrapment.  The  final  formula  for  Pcapt 
assumes  the  form  of 


9.3.  Results  for  Mercury 

For  a  given  quality  function,  the  integral  in  the  Eq.  (38  )  can  be 
computed  numerically  using  the  formulae  (32)  and  (28).  The 
expression  (38)  gives  a  fast  way  of  semi-analytical  estimation  of 
capture  probabilities,  and  it  helps  us  to  avoid  computer-heavy  sim¬ 
ulations.  The  results  of  the  semi-analytical  estimation  of  capture 
probability  as  a  function  of  eccentricity  are  presented  in  Fig.  9  for 
the  set  of  default  parameters  from  Table  2.  The  left  panel  of  the 
graph  shows  the  predicted  probabilities  for  a  solid  uniform  planet 
without  a  core.  The  right  panel  shows,  for  the  same  resonances,  our 
calculations  for  a  realistic  model  of  the  present-day  Mercury  with 
Cm/(Cm  +  Cc)  =  0.5  and  k/Cm  =  1CT5.  This  value  of  the  coefficient 
of  friction  corresponds  to  a  kinematic  viscosity  of  0.008  cm2  s_1, 
which  is  close  to  the  value  for  water  at  30  °C.  We  conclude  that 
even  a  moderate  degree  of  friction  in  the  core-mantle  boundary 
boosts  the  probabilities  of  capture  into  the  3:2  and  higher  reso¬ 
nances.  Capture  into  the  2:1  resonance,  for  example,  becomes  cer¬ 
tain  at  the  current  value  of  eccentricity.  Furthermore,  for  e  =  0.206, 
employment  of  such  a  model  for  Mercury  makes  it  more  likely  for 
the  planet  to  be  trapped  in  the  5:2  resonance  (probability  0.58) 
than  to  traverse  it. 

To  verify  the  semi-analytical  estimate,  we  performed  40  inte¬ 
grations  of  the  Eqs.  (26),  starting  with  a  uniformly  distributed 
7m(0)  and  ym(0)  =  2.009n,  for  e  =  0.1  and  Cm/(Cm  +  Cc)  =  0.5. 
The  equilibrium  spin  rate  of  the  core  for  the  2:1  resonance  is 
yc(0)  =  2.02 n.  The  time  span  of  each  run  was  3.3  x  104  yr,  and 
the  maximal  step  was  2.3  x  1CT3  yr.  The  40  numerical  histories 
rendered  us  19  captures  and  21  passages.  The  probability  of 
capture  estimated  semi-analytically  is  0.48  (see  the  right  panel  in 
Fig.  9).  Thus,  our  semi-analytical  method  and  the  spot-check 
simulation  are  in  excellent  agreement. 


10.  Revisiting  the  Option  3:  Mercury  initially  retrograde 

The  scenario  developed  by  Wieczorek  et  al.  (2012)  comprised 
two  steps: 


P 


capt 


2 

7iVCm 

(Cm+Cc)  (  W(ym)dym+2V2kw0 


(38) 


The  “minus”  sign  in  the  denominator  accounts  for  V  being  negative 
for  supersynchronous  resonances.  By  the  expression  (38),  capture 
probabilities  can  be  estimated  semi-analytically  for  any  resonance 
above  1:1  and  for  any  combination  of  Cm/Cc,  e  and  k. 


Step  1 .  Mercury  begins  with  a  retrograde  spin  and  gets  acceler¬ 
ated  by  tides  (possibly,  assisted  by  the  core-mantle  friction), 
until  trapped  into  the  synchronous  resonance. 

Step  2.  An  impact  disrupts  the  synchronous  spin,  and  Mercury 
shifts  to  the  3:2  resonance. 

Our  tidal  model  prohibits  pseudosynchronous  equilibria.  So,  if 
an  impact  is  energetic  enough  to  disrupt  the  synchronous  spin, 


eccentricity 


Fig.  9.  Probabilities  of  capture  in  the  3:2,  2:1,  5:2  states,  with  the  parameters  as  in  Table  2.  Left:  a  uniform  solid  composition.  Right:  a  molten  core  with  Cm/(Cm  +  Cc)  =  0.5 
and  k/Cm  =  10  5. 
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Fig.  10.  Sample  CDF  (broken  line)  and  expected  population  CDF  (smooth  line)  of  angular  distances  of  large  and  confidently  detected  impact  craters  on  the  surface  of  Mercury 
from  (a)  the  presumably  subsolar  direction  (0°,0°)  and  (b)  East  direction  (90°, 0°). 


but  not  to  reach  the  3:2  resonance,  the  tidal  torque  will  drive  Mer¬ 
cury  into  a  spin-orbit  resonance.  If  the  torque  is  accelerating,  the 
planet  will  be  trapped  into  the  3:2  state.  We  however  explain  in 
the  end  of  this  section  that  the  most  probable  outcome  is  a  nega¬ 
tive  tidal  torque  which  will  return  Mercury  back  into  synchronism. 

To  boost  the  rotation  rate  above  3n/2,  the  impactor  must  have 
created  a  crater  with  a  diameter  of  at  least  600  km.  In  such  a  case, 
we  are  back  to  either  to  the  Option  1  (a  prograde  rigid  planet,  see 
Section  8)  or  to  the  Option  2  (an  initially  prograde  Mercury  with  a 
core,  see  Section  9),  dependent  on  whether  the  core-mantle  fric¬ 
tion  is  present.  We  have  seen  that  in  such  a  case  entrapment  into 
the  3:2  resonance  is  the  most  realistic  outcome.  So,  we  here  just 
need  to  reinvestigate  the  acceleration  of  an  initially  retrograde 
Mercury  with  our  tidal  model,  and  with  or  without  core-mantle 
friction. 

10.1.  Statistical  hypothesis  testing  on  crater  distribution 

To  support  the  hypothesis  of  a  previous  stay  of  Mercury  in  the 
1:1  resonance,  Wieczorek  et  al.  (2012)  looked  for  irregularities  in 
the  distribution  of  large  craters.  Those  counts  were  based  on  lim¬ 
ited  data  obtained  by  Mariner  10  and  two  flybys  of  MESSENGER. 
In  the  light  of  more  accurate  and  complete  data  currently  available, 
a  re-examination  of  the  statistical  confidence  of  the  said  irregular¬ 
ities  is  warranted.  Here  we  use  the  data  published  by  Fassett  et  al. 
(2012,  Table  2)  for  certain  and  verified  craters  larger  than  300  km 
in  diameter.  We  use  the  classical  and  robust  Kolmogorov-Smirnov 
test  of  statistical  hypotheses,  to  estimate  the  confidence  level  of 
any  irregularities  in  the  observed  distribution.  Therefore,  ours  is 
a  completely  independent  verification  of  the  result  published  by 
Wieczorek  et  al.  (2012). 

If  a  distribution  of  points  on  a  unit  sphere  is  absolutely  random 
and  uniform,  the  angular  distance  d  of  each  point  from  a  fixed 
direction  should  follow  a  probability  density  function  (PDF)  that 
is  proportional  to  sin(^)  on  the  range  [0,7i].  Any  significant  non¬ 
randomness  in  the  sample  distribution  will  result  in  a  deviation 
of  the  sample  frequency  of  S  from  this  PDF.  The  corresponding 
cumulative  distribution  function  (CDF)  is  1/2(1  -  cos  5).  The  one- 
sample  Kolmogorov-Smirnov  (KS)  test  determines  the  largest 
deviation  of  the  observed  sample  CDF  from  the  expected  CDF  and 
estimates  the  probability  of  such  deviation  to  happen  by  chance. 

As  explained  in  Wieczorek  et  al.  (2012),  the  most  prominent 
consequence  of  an  extended  period  of  synchronised  rotation  is  a 
bipolar  excess  of  large  craters  on  the  caps  facing  and  opposing 
the  Sun,  i.e.,  close  to  the  directions  (0,0)  and  (180,0)  in  Mercurian 


longitude  and  latitude.  These  directions  are  aligned  with  the  axis  of 
the  smallest  moment  of  inertia  A.  Mathematically,  the  distribution 
of  angular  distances  <5  from  this  chosen  direction  should  be  signif¬ 
icantly  flatter  than  the  expected  PDF.  The  tendency  of  a  synchron¬ 
ised  planet  to  be  bombarded  mostly  from  the  Sun’s  direction  and 
the  opposite  direction  is  caused  by  the  large  fraction  of  asteroids 
on  high-eccentricity  orbits  in  the  population  of  impactors  that 
can  collide  with  the  planet  with  an  impact  velocity  sufficiently 
high  to  leave  behind  a  large  crater. 

The  null  hypothesis  to  test  is  that  the  observed  sample  distribu¬ 
tion  of  S  from  the  fixed  direction  (0,0)  follows  the  general  CDF. 
This  is  done  by  a  one-sample  KS  test  using  standard  statistical  rou¬ 
tines.  The  output  of  the  test,  called  p-value,  quantifies  the  probabil¬ 
ity  of  the  null  hypothesis  to  be  rejected  while  being  true,  given  the 
observed  sample.  In  astronomy  it  is  also  called  false  alarm  proba¬ 
bility  (FAP).  A  small  p-value  indicates  that  it  is  unlikely  that  a  given 
observed  sample  is  drawn  from  a  random,  uniform  distribution.  In 
some  fields,  such  as  exoplanet  detection,  a  high-confidence  result 
is  required  to  have  a  p-value  as  small  as  0.01.  Inversely,  a  large 
p-value  indicates  that  there  is  no  reason  to  reject  the  null  hypoth¬ 
esis.  The  p-value  obtained  from  the  KS  test  for  the  presumably  sub¬ 
solar  point  (0°,  0°)  is  0.946,  so  the  null  hypothesis  is  passable.  There 
is  no  significant  excess  of  crater  density  in  the  direction  aligned 
with  the  minimal-inertia  axis,  as  the  hypothesis  of  synchronous 
spin  predicts.  Fig.  10a  shows  the  observed  sample  CDF  superim¬ 
posed  with  the  null-hypothesis  CDF  for  a  perfectly  random-uni¬ 
form  distribution.  Indeed,  the  two  curves  are  quite  close. 

The  negative  result  for  the  main  prediction  of  the  synchronisa¬ 
tion  surmise  does  not  imply  that  the  observed  distribution  is  ran¬ 
dom  in  any  chosen  direction.  To  test  if  there  are  other  directions  in 
which  the  observed  distribution  may  be  non-random,  we  per¬ 
formed  a  KS  test  for  each  node  of  a  grid  of  1146  points  evenly  dis¬ 
tributed  on  the  sphere.  The  smallest  p-value,  0.006,  was  found  for 
the  coordinates  (11 1.0°,  8. 2°).  Such  a  small  value  is  not  necessarily 
a  good  reason  to  suspect  that  the  sample  is  non-random.  1  5  What  is 
alarming  is  the  closeness  of  this  direction  to  the  due  East  point, 
(90°,  0°).  Indeed,  the  secondary,  and  less  important,  prediction  of 
the  synchronisation  surmise  is  a  dipole-type  asymmetry  in  the 
East-West  distribution  of  craters.  It  is  the  “runner  in  the  rain"  effect, 
when  a  person  gets  soaked  faster  on  the  front  side  than  on  the  back. 
One  can  rationalise  that  the  due  West  direction  is  an  arbitrarily 


15  A  self-consistent  hypothesis  test  requires  the  direction  to  be  chosen  arbitrarily. 
We,  however,  are  dealing  here  with  the  direction  minimising  the  p-value. 
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xM  -  5  yr 


Fig.  11.  Final  states  for  an  initially  retrograde  Mercury  and  rM  =  =  5  yr.  For  a 

longer  Maxwell  time,  we  get  only  synchronous  resonances.  We  can  see  that  the 
trapping  into  the  1:2  is  certain  for  an  eccentricity  e  above  rjO.12. 


chosen  point  of  reference,  not  a  priori  related  to  the  sample,  so  the  KS 
test  is  valid.  The  KS  p-value  in  this  case  is  0.02.  The  corresponding 
pair  of  a  sample  and  expected  CDF  are  displayed  in  Fig.  10b.  The  con¬ 
clusion  is  clear:  there  are  considerably  more  of  large  and  confidently 
detected  craters  on  the  West  side  than  on  the  East. 

Our  statistical  hypothesis  tests  indicate  that  the  primary  predic¬ 
tion  of  the  synchronous  rotation  theory  is  not  supported  by 
observation.  The  secondary  predicted  effect,  i.e.,  the  East-West 
asymmetry  in  the  counts,  is  present  at  a  non-negligible  confidence 
level  of  98%. 

30.2.  A  rigid  and  homogeneous  Mercury 

As  in  Section  8,  we  ran  1000  trajectories  with  the  Maxwell 
times  of  respectively  500.15  and  5  years;  but  now  starting  from 
an  initially  retrograde  spin  corresponding  to  a  period  of  20  days. 
In  this  case,  the  planet  has  to  traverse  spin-orbit  resonances  in 
the  reverse  order  (i.e.,  —2:1,  -3:2,  —1:1,  —1:2,  1:2,  and  1:1),  as 
can  be  understood  from  the  expression  (10)  for  the  triaxiality- 
caused  torque.  The  magnitude  of  the  restoring  component  of  the 
triaxial  torque  is  very  important  for  the  capture  process.  We  recall 
that  the  eccentricity-dependent  coefficients  of  the  restoring 
components,  G2o<j'(e),  scale  in  the  leading  order  as  e4/24, 
e3/48, -e/2. 1,  and  7e/2,  for  q'  =  (-4,  -3, -1, 0, 1),  respectively. 
Obviously,  the  first  two  resonances,  —1:1  and  —1:2,  are  intrinsi¬ 
cally  very  weak  for  moderate  eccentricities,  and  can  in  practice 
be  ignored.  It  is  quite  unlikely  that  planets  can  ever  be  captured 
into  retrograde  spin-orbit  resonances.  The  3:2  resonance  becomes 
stronger  than  the  1 : 1  resonance  at  an  eccentricity  as  small  as  ~0.3. 
The  prograde  (but  subsynchronous)  1:2  resonance  can,  however, 
become  significant,  depending  on  the  eccentricity.  At  first  glance, 
this  resonance  should  correspond  to  an  unstable  equilibrium, 
because  the  corresponding  coefficient  C20(-i)(e)  is  negative  and 
the  torque  becomes  diverging  rather  than  restoring.  However, 
our  simulations  have  shown  that,  under  certain  circumstances, 
capture  into  this  resonance  is  probable  or  even  unavoidable.  As 
mentioned  by  Goldreich  and  Peale  (1966),  this  resonance  is  stable 
if  the  planet  is  captured  sidewise,  i.e.,  if  its  shortest  axis  of  inertia  is 
at  7t/2  to  the  Sun  direction  at  periastron. 

For  tm  =  ta  =  500  yr  and  15  years,  all  the  simulated  trajectories 
were  eventually  trapped  into  the  synchronous  resonance, 
successfully  crossing  the  1:2  resonance.  However,  the  results  for 
tm  =  zA  =  5  yr  were  quite  different,  with  878  trajectories  ending 
up  in  the  1:2  resonance  (see  Fig.  11).  Thus,  if  the  early  Mercury 
was  warm  and  had  a  low  viscosity,  it  could  hardly  avoid 


entrapment  into  this  sub-synchronous  resonance.  The  simulations 
also  showed  that,  once  captured,  the  planet  would  stay  in  this  res¬ 
onance  indefinitely  long,  despite  significant  variations  in  eccentric¬ 
ity.  Only  one  history  spontaneously  left  the  1:2  resonance.  (It 
ended  up  in  the  1:1  one).  The  conclusion  is  that  a  warm  Mercury 
with  an  initial  retrograde  spin  would  have  an  extremely  slim 
chance  of  getting  to  synchronism.  En  route  thereto,  it  would  almost 
certainly  have  been  trapped  into  the  1:2  subsynchronous  state. 
Only  external  agents,  such  as  giant  impacts,  could  drive  it  through 
to  the  current  3:2  resonance.  But  the  likely  capture  into  the  1:2 
resonance  becomes  another  problem  for  the  retrograde  spin 
hypothesis,  because  at  least  two  fortuitous  giant  impacts  are  now 
needed  to  account  for  the  current  state. 


3  0.3.  A  segregated  Mercury  with  a  fluid  core 


We  have  not  conducted  large-scale  numerical  simulations  of  a 
differentiated-Mercury  spin,  because  the  analysis  presented  in  Sec¬ 
tion  9  is  fully  applicable  in  this  case  too.  Undoubtedly,  a  segregated 
planet  with  a  massive,  decoupled  liquid  core  would  be  even  more 
easily  captured  into  the  subsynchronous  1:2  resonance  than  its 
homogeneous  analogue. 

Examining  the  hypothesis  that  Mercury  was  entrapped  into 
synchronism  for  a  continuous  length  of  time,  and  that  bombard¬ 
ment  by  heavy  asteroids  drove  it  out  to  the  present  3:2  spin  state, 
one  has  to  enquire  whether  at  that  time  the  planet  already  had  a 
liquid  core.  Scenarios  involving  a  synchronous  state  were  exam¬ 
ined  numerically  by  Correia  and  Laskar  (2012),  in  the  framework 
of  the  CTL  tidal  model.  Here  we  offer  analysis  more  relevant  for 
our  Maxwell-Andrade  tidal  model. 

Bringing  a  synchronous  planet  into  a  3:2  spin-orbit  resonance 
comprises  two  steps.  First,  the  planet  should  be  kicked  out  of  syn¬ 
chronism  by  a  mega-impact.  For  small  and  moderate  eccentricities 
(e  <  0.3),  the  1:1  resonance  is  the  most  stable  equilibrium.  It  is  the 
deepest  well  of  potential  energy;  to  leave  it,  a  formidable  external 
energy  deposition  is  needed.  The  direction  of  impact  should  be 
aligned  with  the  sidereal  rotation  velocity  (lest  Mercury  finds  itself 
in  a  sub-synchronous  state  whence  it  is  destined  to  fall  into  either 
1:2  or  1:1  resonance  again).  Second,  Mercury’s  spin  should  evolve 
toward  the  3:2  resonance  after  the  impact.  The  state  of  pseudosyn¬ 
chronism  at  spin  rates  between  the  1:1  and  3:2  resonance  provides 
a  convenient  and  stable  transit  point  in  the  CTL,  but  is  unstable  for 
planets  of  terrestrial  composition. 

There  remain  two  possible  scenarios  consistent  with  the  impact 
hypothesis.  First,  the  impact  could  be  energetic  enough  to  drive 
Mercury’s  spin  rate  to  values  above  the  3:2  resonance,  0  >  3/2 n. 
We  can  estimate  the  characteristics  of  the  required  impactor.  If  u 
is  the  component  of  the  relative  impact  velocity  in  the  equatorial 
plane  of  Mercury  (supposed  to  be  close  to  the  orbital  plane,  due 
to  a  small  pre-impact  obliquity),  the  upper-bound  change  of  the 
spin  rate  from  an  approximation  ignoring  possible  ejecta  and  ther¬ 
mal  dissipation  is 


Am  sC 


m,uR 

C 


(39) 


m,  being  the  impactor’s  mass,  while  C  and  R  being  the  principal 
moment  of  inertia  and  the  radius  of  Mercury.  In  our  approximation, 
the  upper  bound  is  achieved  when  the  collision  is  grazing  and  on 
the  equator.  On  the  other  hand,  the  change  of  spin  rate  should  be 
greater  than  the  gap  1  /2 n  separating  the  two  resonances.  The 
upper-bound  condition  becomes 

(40) 


For  the  present-day  parameters,  the  momentum  of  the  impactor 
should  have  exceeded  1.13  x  1023  kg  m  s_1,  if  Mercury  were  a  uni- 
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formly  rigid  body.  According  to  Le  Feuvre  and  Wieczorek  (201 1 ),  the 
average  impact  velocity  is  42.5  km  s  The  average  projection 
component  u  is  statistically  smaller,  2/7t  ■  42.5  =  27.1  km  s_1.  The 
minimum  mass  of  the  impactor  capable  of  driving  a  uniform 
Mercury  analogue  all  the  way  above  the  3:2  resonance  is  then 
4.2  x  1018  kg.  Using  an  assumed  median  density  of  asteroids  of 
2500  kg  nr3,  and  an  estimated  relation  between  the  basin  and 
impactor  diameters,  we  conclude  that  only  the  largest  two  of  the 
certain  and  probable  craters  from  Fassett  et  al.  (2012)  -  namely, 
Caloris  and  b30  -  could  impart  that  much  momentum  under 
favourable  circumstances.  The  presence  of  a  large  liquid  core  decou¬ 
pled  from  the  mantle  lowers  the  threshold  momentum  and  the 
mass  of  impactor  by  roughly  one-half,  as  only  the  moment  of  the 
mantle’s  inertia,  Cm,  enters  the  Eq.  (40).  Consequently,  any  impactor 
leaving  a  basin  larger  than  ~640  km  could  potentially  cause  the 
required  spin-up.  Still,  there  are  only  a  dozen  or  so  such  craters 
on  Mercury. 

In  the  alternative  scenario,  the  planet  only  has  to  be  driven  out 
of  the  1:1  resonance,  the  subsequent  spin-up  and  capture  into  the 
3:2  resonance  being  assured  by  tides.  The  required  impactor’s 
momentum  is  less  stringent  then.  The  lower-bound  trajectory  of 
departure  from  the  1:1  resonance  is  obviously  the  same  as  the 
zero-energy  trajectory  of  capture,  i.e.,  the  separatrix  from  the 
Eq.  (28).  The  corresponding  lower-bound  change  of  spin  rate  is 

(41) 

ignoring  the  small  forced  libration  which  is  slightly  dependent  on 
the  orbital  eccentricity  in  this  case.  We  obtain  the  lower-bound 
Ao)mi„  of  2.0  x  10  s,  1.9  x  HT8,  and  1.8  x  1CT8  Hz,  for  e  =  0.1,  0.2 
and  0.3,  respectively.  The  lower-bound  impactor’s  momentum  is 
then 
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and  the  corresponding  minimum  mass  ranges  between  1.8  x  1017 
and  2.0  x  1017  kg.  The  corresponding  minimal  basin  diameter  is 
approximately  510  km.  Thus,  half  of  the  impactors  considered  in 
Section  10.1  could  do  the  work  of  driving  the  planet  out  of  synchro¬ 
nism  under  favourable  conditions.  Again,  the  required  impactor’s 
size  becomes  even  smaller  for  a  planet  with  a  decoupled  core.  The 
lower  bound  of  basin  diameter  goes  down  to  ~430  km,  which  adds 
several  more  potentially  capable  impact  event,  for  which  we  have 
observational  evidence. 

The  main  obstacle  for  this  scenario  to  be  accepted  as  plausible, 
with  or  without  a  liquid  core,  is  the  necessity  to  further  accelerate 
the  spin  from  the  post-impact  value  to  the  3:2  commensurability. 
Within  our  realistic  tidal  model,  we  do  not  have  the  help  from  a 
stable  intermediate  equilibrium  of  pseudosynchronous  rotation. 
Thus,  the  required  spin-up  should  be  provided  by  a  steady  positive 
tidal  torque.  This  imposes  a  requirement  on  the  value  of  eccentric¬ 
ity  during  this  time  span.  For  each  non-resonant  value  of  0,  there  is 
a  unique  equilibrium  value  of  eccentricity  that  separates  the 
domains  of  negative  and  positive  secular  tidal  torque.  The 
equilibrium  eccentricity  for  a  uniform  Mercury  analogue  was  com¬ 
puted  in  Makarov  (2012,  Fig.  3).  The  area  of  interest  here  is  the 
span  of  0  between  In  and  1 ,5n  where  the  eccentricity  is  descending 
from  approximately  0.3  to  approximately  0.2.  The  phase  space  tra¬ 
jectory  of  the  planet  depends  on  the  initial  post-impact  spin  rate.  If 
the  impact  was  just  enough  to  drive  the  planet  out  of  the  1:1  res¬ 
onance,  the  initial  spin  rate  would  end  up  right  above  1  n,  in  which 
case  the  eccentricity  would  have  to  be  consistently  above  0.3,  to 
boost  the  spin  further.  Otherwise,  Mercury  quickly  returns  to  the 
1:1  resonance.  A  more  powerful  impact  could  set  the  spin  rate  at 
a  higher  value,  which  eases  up  the  limit  on  eccentricity  needed 


for  a  further  prograde  spin-up.  Statistical  analysis  can  be  per¬ 
formed  by  taking  into  account  the  observed  distribution  of  impac¬ 
tors  and  the  Monte-Carlo  simulations  of  eccentricity  evolution,  but 
such  an  endeavour  is  outside  the  scope  of  this  paper. 


11.  Discussion 

A  large  variety  of  scenarios  of  Mercury’s  spin-orbit  evolution 
has  been  considered  in  the  literature,  striving  to  explain  its  present 
3:2  resonance.  Differentiating  the  proposed  theories  by  the  degree 
of  complexity  and  physical  validity,  the  options  include:  (1)  simpli¬ 
fied  linear  tidal  models  versus  physics-based  rheologies  resulting 
in  realistic  frequency-dependencies  of  k2/Q_\  (2)  constant-in-time 
versus  chaotically  variable  eccentricity;  (3)  initially  prograde 
versus  retrograde  spin;  (4)  the  possibility  of  powerful  impacts  fol¬ 
lowed  by  instantaneous  resets  of  the  spin  rate;  (5)  the  presence  or 
absence  of  a  massive  molten  core  at  the  epoch  of  capture  into  the 
3:2  resonance;  (6)  the  tidal  properties,  such  as  the  average  viscos¬ 
ity,  of  Mercury’s  mantle  at  the  time  of  capture.  Thus  a  multitude  of 
possible  paths  emerges,  and  the  results  happen  to  strongly  depend 
upon  the  choice  made  on  each  of  the  said  six  items.  Not  all  of  these 
paths,  however,  render  the  current  3:2  resonance  as  the  likeliest 
end-state. 

Reliable  conclusions  on  the  history  of  Mercury’s  spin  must  be 
derived  from  the  actual  behaviour  of  minerals  -  a  viscoelastic  Max¬ 
well  response  at  low  frequencies  and  the  inelastic  Andrade  creep  at 
higher  frequencies.  This  complicates  both  numerical  simulations 
and  direct  analysis  of  capture  probabilities.  Another  complication 
is  the  large,  unpredictable  variation  of  eccentricity  over  aeons. 
The  onus  is  then  on  massive  Monte-Carlo  simulations,  though  we 
also  found  a  limited  analytical  estimate  to  be  valuable  as  a  sanity 
check.  In  other  regards,  a  wide  scope  of  initial  conditions  and 
assumptions  has  been  considered  in  this  paper,  avoiding  a  priori 
preferences.  We  explored  both  the  prograde  and  retrograde  initial 
spin;  considered  a  planet  with  or  without  a  liquid  core  at  the  early 
stage  of  its  evolution;  and  studied  both  a  partially  molten  warm 
mantle  and  a  colder  mantle  similar  to  that  of  the  present-day 
Earth. 

The  main  outcome  of  our  study  is  that  the  3:2  spin-orbit  reso¬ 
nance  is  the  likeliest  end-state,  provided  the  following  set  of  con¬ 
ditions  is  satisfied  by  Mercury:  (a)  it  did  not  have  a  massive 
decoupled  core  prior  to  being  captured  into  this  resonance;  (b) 
its  mantle  was  relatively  cold  and  had  a  sufficiently  high  viscosity, 
with  a  Maxwell  time  being  of  the  order  of  hundreds  years;  (c)  the 
initial  spin  was  prograde,  as  would  be  expected  for  most  of  the 
Solar  system  planets.  An  initially  prograde  undifferentiated  cold 
Mercury  gets  trapped  into  the  resonance  within  10-20  Myr,  while 
formation  of  a  core  requires  up  to  a  billion  of  years. 

Our  analysis  thus  invalidates  some  previously  discussed 
scenarios.  We  do  not  find  it  possible  for  a  telluric  planet  locked 
in  a  spin-orbit  resonance  to  leave  it  spontaneously  through  a 
change  in  eccentricity.  In  our  simulations,  resonant  states  were 
almost  always  long-term  and  stable.  Without  an  external  agent, 
such  as  a  heavy  impact,  Mercury  is  trapped  forever.  This  fact,  com¬ 
bined  with  the  short  spin-down  times  found  in  this  paper,  implies 
that  Mercury  was  captured  into  the  3:2  state  billions  of  years  ago, 
and  stayed  there  ever  since. 

Alternatives  to  this  scenario  meet  substantial  difficulties.  If  the 
initial  spin  were  fast  and  prograde,  but  the  core  had  already  been 
segregated  or/and  if  the  mantle  was  warm  by  the  moment  of  cap¬ 
ture,  then  Mercury  would  very  likely  end  up  in  a  higher  resonance, 
e.g.,  2:1.  If  the  initial  spin  were  retrograde,  then  the  hitherto 
unheeded  sub-synchronous  resonance  1:2  becomes  hard  to  avoid, 
especially  for  a  differentiated  or/and  warm  planet.  So  two  hits 
would  be  needed  -  to  kick  the  planet  from  the  likely  trap  of  1:2, 
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and  then  to  kick  it  out  from  the  certain  trap  1:1.  Hence,  both  for  an 
initially  prograde  and  an  initially  retrograde  planet,  early  differen¬ 
tiation  becomes  a  formidable  obstacle  on  its  way  to  the  3:2  spin 
state. 

The  published  studies  of  Mercury’s  differentiation  are  centred 
around  the  fact  that,  being  anomalously  dense,  the  planet  should 
be  poor  in  silicates.  Early  volatilisation  is  one  of  the  attractive  pos¬ 
sibilities.  The  composition  can  be  explained  if  the  removal  process 
in  Mercury’s  zone  of  the  solar  nebula  were  only  slightly  more  effec¬ 
tive  for  silicates  than  for  iron.  The  fractionation  required  to  pro¬ 
duce  an  iron-rich  planet  is  achieved  through  a  combination  of 
gravitational  and  drag  forces  acting  in  the  nebula  during  the 
formation  of  proto-Mercury  (Weidenschilling,  1978).  As  an  alter¬ 
native,  a  slow-volatilisation  model  was  suggested,  in  which  the 
solar  wind,  after  the  formation  of  Mercury,  could  partially  volatil¬ 
ise  the  mantle  (Cameron,  1985;  Fegley  and  Cameron,  1987).  Both 
these  models  imply  gradual  differentiation  of  the  planet.  As 
explained  in  Appendix  A,  the  process  of  core  formation  required 
hundreds  of  millions  to  a  billion  of  years.  Simulations  without 
core-mantle  friction  furnishes  entrapment  within  less  than 
20  Myr,  an  age  indeed  much  younger  than  that  of  differentiation. 

A  very  different  model  of  rapid  differentiation  relies  on  an 
impact  energetic  enough  to  melt  down  the  originally  chondritic 
proto-Mercury  (Benz  et  al.,  2007).  The  event  would  have  resulted 
in  a  very  warm  tumbling  Mercury  wherein  the  core-mantle 
boundary  would  have  formed  very  quickly.  An  impact,  which  could 
have  taken  place  when  Mercury’s  age  was  about  20  Myr  or  larger, 
is  consistent  with  our  core-less  computation:  the  capture  happens 
first,  the  impact  and  the  subsequent  differentiation  comes  later.  An 
earlier  impact  poses  a  nontrivial  problem  for  our  main  scenario. 
The  hypothesis  of  an  early  impact  becomes  questionable  also  in 
the  light  of  the  recent  findings  by  MESSENGER.  As  was  pointed 
out  by  Peplowski  et  al.  (2011),  the  abundance  of  potassium  is 
inconsistent  with  those  physical  models  of  Mercury  formation,  which 
require  extreme  heating  of  the  planet.  This  observation  prompted 
Wurm  et  al.  (2013)  to  explore  the  action  of  photophoretic  forces 
upon  irradiated  solids.  The  investigation  has  demonstrated  that  sil¬ 
icates  are  preferentially  pushed  into  an  optically  thick  disc.  The 
subsequent  planetesimal  production  at  the  outward-moving  edge 
renders  metal-rich  platenesimals  close  to  the  star  and  metal- 
depleted  ones  farther  out  in  the  nebula.  After  this  very  early  sepa¬ 
ration  of  metals  and  silicates,  the  further  development  of  Mercury 
goes  gradually,  with  a  slow  emergence  of  a  core-mantle  boundary. 

A  massive  impact  seems  to  provide  a  magic  resolution  to  many 
of  these  difficulties.  Wieczorek  et  al.  (2012)  suggest  that  the  distri¬ 
bution  of  craters  indicates  that  Mercury  spent  some  time  in  the 
synchronous  state  whence  it  was  driven  by  an  impact,  eventually 
reaching  the  current  3:2  resonance.111  We  reconsidered  the  claim 
of  inhomogeneous  distribution  of  impact  craters  on  the  surface  of 
Mercury,  using  more  up-to-date  and  complete  data  and  a  robust 
method  of  statistical  hypothesis  testing.  No  evidence  is  found  for  a 
preferential  location  of  craters  on  the  opposite  caps  in  the  sub-solar 
and  anti-solar  directions,  which  is  the  principal  prediction  from  the 
synchronous-state  hypothesis.  Still,  the  secondary  prediction  of  an 
East-West  asymmetry  in  the  distribution  of  craters  is  confirmed. 

A  consistent  tidal  theory  based  on  rheology  of  minerals  is,  by 
itself,  insufficient  to  disprove  the  scenario  by  Wieczorek  et  al. 
(2012).  The  scenario,  however,  comprises  a  number  of  propositions 


16  Wieczorek  et  al.  (2012)  hypothesise  that  Mercury  could  have  evolved  into 
synchronism  from  an  originally  retrograde  spin  state.  The  probability  of  this  step  is 
reduced  considerably  by  the  likely  entrapment  into  the  1:2  resonance,  escape 
wherefrom  would  require  an  extra  collision.  At  the  same  time,  Brasser  et  al.  (2013) 
argue  that  the  primordial  orbit  of  Mercury  could  have  been  more  circular  than 
previously  thought.  This  finding,  if  correct,  may  give  the  second  breath  to  the  scenario 
with  synchronism.  If  the  eccentricity  keeps  a  veiy  low  value  for  several  millions  of 
years,  a  swift  despinning  straight  into  synchronism  is  not  impossible. 


which  are  not  impossible  when  considered  separately,  but  which 
together  combine  into  an  unlikely  course  of  events.  Within  this 
scenario,  an  initially  retrograde  spin  slowed  down  and  became 
prograde.  Then  Mercury  would  likely  get  into  the  first  sub- 
synchronous  resonance  1:2,  an  event  almost  inevitable  for  inviscid 
planets  with  a  significant  liquid  core.  Mercury  would  hardly  be 
able  to  leave  this  state  -  hence  the  need  for  a  heavy  impact  in 
the  right  direction.  The  probability  that  a  single  impact  is  directed 
in  the  prograde  sense  is  1  /2.  A  counter  strike  could  reset  the  spin 
back  to  smaller  values,  whereafter  Mercury  would  fall  in  the  1:2 
state  again.  Once  externally  accelerated  toward  the  1:1  resonance, 
Mercury  would  reach  it  in  a  matter  of  millions  of  years  and  get 
inevitably  captured  there,  irrespective  of  its  eccentricity.  Only  most 
energetic  impacts  (again,  correctly  directed  and  almost  tangential 
to  the  surface  in  the  equatorial  zone)  can  drive  Mercury  all  the 
way  to  spin  rates  above  the  3:2  resonance.  There  are  few  craters 
that  could  be  left  by  such  massive  impactors,  so  this  option  appears 
to  be  unlikely.  It  is  much  easier  to  jolt  Mercury  out  of  the  1:1  state 
to  spin  rates  just  outside  the  resonance,  especially  if  it  already  pos¬ 
sesses  a  large  liquid  core.  Still,  if  the  eccentricity  at  the  moment  is 
smaller  than  0.3,  Mercury  would  quickly  descend  back  into  syn¬ 
chronism.  Again,  any  counter-directed  strike  would  only  enable  a 
short  excursion  outside  the  main  resonance.  For  Mercury  to  get 
to  the  3:2  state,  being  driven  by  the  tides  only,  the  eccentricity 
should  be  higher  than  0.3-0.2  for  the  whole  duration  of  spin-up. 

The  scenario  of  initially  prograde  Mercury  with  a  core  meets  its 
own  set  of  difficulties  in  our  analysis.  Given  the  short  despinning 
time,  the  emergence  of  a  core  prior  to  capture  is  unlikely.  (Tidal 
damping  in  Mercury  is  too  weak  to  account  for  a  high  degree  of  seg¬ 
regation  Makarov  and  Efroimsky  (2014).)  Besides,  a  decoupled  core 
makes  capture  into  the  2:1  resonance  almost  certain,  and,  as  our 
simulations  show,  permanent.  Hence  we  have  to  resort  to  the  agency 
of  an  external  impact  to  explain  the  current  3:2  resonance.  This  is 
not  entirely  impossible,  simply  because  it  is  much  easier  to  drive 
Mercury  down  from  the  2:1  resonance  than  up  from  the  1:1  reso¬ 
nance,  with  small  and  moderate  eccentricities.  From  the  energy-bal¬ 
ance  viewpoint,  it  is  more  economic  to  descend  down  a  bumpy  slope 
than  to  be  hoisted  out  of  a  hole  at  the  bottom.  Altogether,  though, 
this  scenario  requires  too  many  unlikely  assumptions. 

To  draw  to  a  close,  we  advocate  a  capture  of  an  initially 
prograde,  yet  undifferentiated  Mercury.  It  requires  less  than 
20  Myr,  the  3:2  resonance  being  the  most  probable  end-state. 

Acknowledgments 

Benoit  Noyelles  is  an  F.R.S.-FNRS  postdoctoral  research  fellow, 
and  he  thanks  this  institution  for  having  funded  two  of  his  visits 
to  the  USNO.  This  research  has  made  use  of  the  resources  of  the 
Interuniversity  Scientific  Computing  Facility  located  at  the  Univer¬ 
sity  of  Namur  and  supported  by  the  F.R.S.-FNRS  under  Convention 
No  2.4617.07.  Julien  Frouard  acknowledges  support  by  NSF  grant 
AST-1109776. 

On  various  occasions,  the  topics  addressed  in  the  paper  were 
discussed  by  the  authors  with  Sylvio  Ferraz  Mello,  Valery  Lainey, 
Francis  Nimmo,  Stanton  Peale,  James  G.  Williams,  Jack  Wisdom 
and  other  colleagues,  to  all  of  whom  the  authors  are  profoundly 
grateful.  The  authors  are  thankful  to  Sebastiano  Padovan  for  a  con¬ 
sultation  on  the  rescaling  method  employed  in  Padovan  et  al. 
(2014).  The  authors  are  indebted  to  Valery  Lainey  and  Christophe 
Le  Poncin-Lafitte  who  made  this  collaboration  possible  by  putting 
the  authors  in  touch  with  one  another. 

Appendix  A.  How  long  it  takes  Mercury  to  develop  a  liquid  core 

Our  modelling  shows  that,  for  a  realistic  rheology,  capture  hap¬ 
pens  extremely  swiftly:  almost  certainly  within  20  Myr,  more 
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probably  within  10  Myr.  By  that  age,  the  liquid  core  of  Mercury  had 
not  yet  formed.  Indeed,  over  a  time  span  At,  the  temperature 
increases  by 

AT  =  ^-At,  (Al) 

bp 

where  we  neglected  the  diffusion,  and  denoted  with  H  and  Cp  the 
radiogenic  heating  and  specific  heat  capacity,  correspondingly. 
Using  the  present-day  values  H  =  3.5  x  1CT12  W/kg  and 
Cp  =  1200 J/(I<  kg)  borrowed  from  Turcotte  and  Schubert  (2002), 
we  observe  that  by  the  age  of  At  =  20  Myr  the  radiogenic  tempera¬ 
ture  increase  of  Mercury  was  AT  =  1 .84  K  only.  In  the  early  Solar 
System,  the  value  of  H  could  have  been  up  to  four  times  higher, 
but  this  would  be  still  too  low  to  melt  the  mantle  of  Mercury  over 
such  a  short  time.  A  likely  age  of  melting  is  hundreds  of  millions  to  a 
billion  of  years.1 ' 

In  Section  9,  we  explore  trapping  of  a  Mercury  with  a  liquid 
core.  Our  semi-analytical  model  (confirmed  by  limited  numerical 
simulation)  shows  that  even  a  weak  laminar  friction  between  the 
mantle  and  the  core  would  result  in  a  capture  into  the  2:1  or  a 
higher  resonance. 

These  two  circumstances  indicate  that  Mercury  was  homoge¬ 
neous  when  trapped. 

Even  for  a  homogeneous  Mercury,  a  decrease  of  the  Maxwell 
time  from  hundreds  to  dozens  of  years  (which  corresponds,  very 
roughly,  to  a  100  K  temperature  increase)  makes  the  2:1  resonance 
a  more  probable  destination  than  3:2  (see  the  Section  8.2).  A  fur¬ 
ther  decrease  of  the  Maxwell  time  makes  the  3:2  final  resonance 
even  less  probable.  On  all  accounts,  the  warmer  an  early  Mercury 
the  less  chances  it  has  to  end  up  where  it  is  now. 


Appendix  B.  The  separatrix  dividing  circulation  from  libration. 
Derivation  of  the  Eq.  (28) 


In  neglect  of  friction,  the  law  of  motion  near  the  q'  resonance 
can  be  written  as 


Cm7  +  2  (Bm  Am)  M 


-  +  M, 


planet 


-n2G2oq'(e)  sin2y  =  0. 


(Bl) 


As  was  pointed  out  by  Goldreich  (1966),  the  problem  is  identical  to 
a  simple  pendulum.  Multiplying  the  Eq.  (Bl)  by  y  and  integrating 
over  time,  we  get  an  integral  of  motion, 
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An) 


Ms, 


Mstar  Mpianet 


n2G2  oq-(e)  cos  2  y  =  E, 


(B2) 


whose  value  is  defined  by  the  initial  conditions.  There  is  a  value 
E  =  Eb  separating  circulation  from  libration.  To  find  Eb,  we  set 
y  =  7t/2  and  y  =  0  in  the  (B2).  This  furnishes: 


M, 


planet 


-C2oq'(e)n2 


(B3) 


insertion  whereof  into  the  formula  (B2)  renders  us  the  equation  of 
the  separatrix18: 


3 (Bm  -Am) 
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planet 


-G20q'(e) 
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cosy. 


(B4) 


17  The  authors  are  grateful  to  Francis  Nimmo  for  the  help  he  provided  with  this 
estimate. 

18  Our  Eqs.  (Bl)  and  (B2)  and  (B4)  coincide  with  those  in  Goldreich  and  Peale  (1968), 
but  differ  from  similar  equations  in  Makarov  (2012)  and  Makarov  et  al.  (2012).  The  y 
used  in  Goldreich  and  Peale  (1968)  and  in  the  current  paper  differs  by  a  factor  of  1/2 
from  the  y  in  Makarov  (2012)  and  Makarov  et  al.  (2012). 
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